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ABSTRACT 


The  layers  of  the  earth's  crust  act  as  a  filter  with 
respect  to  seismic  energy  arriving  at  a  given  station. 
Consequently  the  motion  recorded  at  the  earth's  surface  de¬ 
pends  not  only  on  the  frequency  content  of  the  exciting 
seismic  energy  and  on  the  response  characteristics  of  the 
recording  instrument,  but  also  on  the  elastic  parameters 
and  thicknesses  of  the  layers.  This  latter  dependence  is 
the  basis  for  a  method  of  investigating  the  structure  of 
the  crust. 

In  order  to  obtain  information  independent  of  the  time 
history  and  spatial  distribution  of  the  source  of  energy 
the  spectrum  of  the  vertical  component  of  motion  Is  divided 
by  the  spectrum  of  the  horizontal  component.  This  ratio 
represents  the  tangent  of  the  apparent  angle  of  emergence 
as  a  function  of  frequency.  It  depends  only  on  the  angle 
of  incidence  of  the  ray  and  the  system  of  layers  below  the 
recording  station.  The  parameters  of  the  crust  may  be  de¬ 
termined  by  comparison  of  theoretical  and  observed  spectra 
of  this  ratio. 

To  facilitate  this  comparison,  a  set  of  master  curves 
was  calculated  using  the  matrix  development  of  Haskell. 
Calculations  of  these  curves  are  in  terms  of  a  dimension¬ 
less  frequency.  This  presentation  allows  the  grouping  of 
the  curves  corresponding  to  different  crustal  models  into 
families  of  curves.  A  set  of  master  curves  of  the  appar¬ 
ent  angle  of  emergence  for  one-layer  models  and  for  dif¬ 
ferent  angles  of  incidence  and  contrasts  of  velocities  be¬ 
tween  the  crust  and  the  mantle  is  presented.  This  set  Is 
complete  In  the  sense  that  any  one-layer  model  may  be 
Interpolated.  A  second  set  for  some  combinations  of  two- 
layer  models  Is  also  presented. 

The  characteristics  of  these  curves  are  discussed  from 
the  point  of  view  of  their  "periodicity"  in  the  frequency 
domain  and  of  their  amplitude  in  order  to  Investigate  the 
Influence  of  the  layer  parameters.  Considerations  either 
of  constructive  interference  or  of  Fourier  analysis  of  a 
pulse  multiply  reflected  within  the  crust  reveal  that  the 
amplitude  of  peaks  and  troughs  in  the  spectrum  depends  on 
the  velocity  contrast  at  the  interfaces  of  the  system.  The 
"periodicity"  or  spacing  of  peaks  and  troughs  depends  on 
the  time  lags  between  the  first  arrival  of  the  direct  P 
wave  and  the  secondary  arrivals  of  the  converted  waves  or 
of  multiply  reflected  and  refracted  waves.  Closely  spaced 
fluctuations  correspond  to  large  time  lags,  and  widely 
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spaced  fluctuations  to  short  time  lags. 

Observations  of  the  spectrum  of  the  apparent  angle  of 
emergence  were  obtained  by  dividing  the  smoothed  spectra 
of  the  vertical  and  horizontal  component  seismograms.  In 
order  to  avoid  the  influence  of  reflections  at  the  crust 
near  the  source  or  of  reflections  from  the  core  of  the 
earth,  the  earthquakes  selected  were  of  intermediate  and 
large  focal  depth  and  were  restricted  to  epicentral  dis¬ 
tances  less  than  55°. 

Application  of  the  method  to  the  long-period  seismo¬ 
grams  of  the  Saint  Louis  University  Network  of  stations 
gives  an  average  P  velocity  in  the  crust  of  6.6  km/sec 
and  a  total  thickness  of  the  crust  of  42  km  for  the  cen¬ 
tral  part  of  the  United  States  under  these  stations. 
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THE  DETERMINATION  OF  CRUSTAL  THICKNESS  FROM  THE 


SPECTRUM  OF  P  WAVES 

1.  Introduction 

For  many  purposes  the  crust  of  the  Earth  and  the 
upper  mantle  may  be  considered  as  a  system  of  horizontal 
layers.  When  this  system  is  excited  by  seismic  energy  its 
frequency  response  to  the  forced  vibration  is  a  function 
of  the  elastic  parameters  and  thicknesses  of  the  layers. 

As  a  result  of  this  behavior  the  seismograms  obtained  at 
the  surface  contain  detailed  information  pertaining  to 
the  crust  and  upper  mantle  of  the  earth  immediately  below 
the  point  of  observation..  Modern  advances  in  the  .heory 
of  elastic  waves  transmitted  by  layered  media  and  computing 
facilities  by  means  of  electronic  computers  make  it  pos¬ 
sible  today  to  extract  this  Information  and  to  evaluate 
to  a  greater  or  lesser  degree  the  parameters  of  the  crust 
from  the  spectrum  of  longitudinal  seismic  waves. 

The  analysis  of  the  body  waves  for  tnis  purpose  Is 
more  conveniently  performed  in  the  frequency  domain.  This 
is  the  result  of  recent  studies  by  Haskell  (1953,  3962), 
Phinney  (1964),  Hannon  (1964b),  Harkrider  (1964)  and  Fuchs 
(1965),  who  have  investigated  the  theoretical  calculation 
of  the  frequency  response  of  layered  systems  to  the  energy 
of  body  waves . 
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In  the  present  investigation  the  effect  of  layered 
media  on  infinite  trains  of  sinusoidal  longitudinal  plane 
waves  incident  obliquely  from  below  is  studied.  The  re¬ 
sponse  of  different  stratified  models  is  systematized  in 
<*uch  a  way  that  families  of  curves  for  different  models  are 
obtained.  By  interpolation  the  response  of  intermediate 
models  is  made  possible.  Since  the  apparent  angle  of  emer¬ 
gence  of  longitudinal  waves  is  Independent  of  the  frequency 
content  of  the  source,  the  effect  of  the  layered  system  on 
this  angle  is  taken  as  the  basis  of  the  analysis. 

1.1  Observations  of  the  effect  of  local  geological  param¬ 
eters  on  longitudinal  seismic  waves. 

Instrumental  observations  of  the  effect  of  the  crustal 
layer  on  the  periods  and  amplitudes  of  seismic  body  waves 
were  first  reported  by  certain  Japanese  authors:  Imamura 
(1929),  Ishimoto  (1931a,  1931b,  1932,  1934),  Nasu  (1931), 
Inouye  (1934)  and  Takahasi  and  Hiramo  (1941).  These  Inves¬ 
tigators  reported  differences  in  the  seismograms  obtained 
for  the  same  earthquake,  at  the  same  eplcentral  distance 
and  azimuth,  and  by  similar  instruments  but  located  on  dif¬ 
ferent  geological  sites.  As  a  possible  explanation  of  the 
dlffere.ee  the  investigators  suggested  the  reflecting  char¬ 
acter  of  the  surface  layers.  At  the  same  time  Suzuki  (1932) 
observed  changes  in  the  apparent  angle  of  emergence  with 
change  in  period  of  the  waves.  His  explanation  of  the 
phenomenon  was  also  in  terms  of  the  reflections  and  refrae- 
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tions  taking  place  within  the  layers  of  the  crust. 

Similar  types  of  observations  were  made  in  Europe  and 
America.  Gutenberg  (1934)  noticed  the  preference  of  cer¬ 
tain  stations  to  detect  waves  of  a  given  period.  The  ampli¬ 
tude  of  the  vibrations  was  also  connected  with  the  geologi¬ 
cal  qualities  of  the  site,  Neumann  (1954)  noticed  the 
different  'ecoense  of  the  '’granitic  type  of  rock"  and  "deep 
sedimentary  rock."  Greater  amplitudes  are  always  connected 
with  the  second  type  of  rock.  The  amplitude  of  the  seismic 
waves  in  relation  to  the  geology  of  the  place  where  the 
instruments  were  located  was  also  studied  by  Gutenberg 
(1957).  He  further  noticed  that  the  long-period  waves  are 
not  so  affected  by  geological  variations  as  are  short-period 
waves.  The  effect  of  the  ground  on  the  amplitudes  of  body 
waves  is  the  main  reason  why  Gutenberg  and  Richter  (1956) 
introduced  the  station  correction  factors  when  the  magnitude 
of  the  earthquakes  was  determined  from  the  amplitude  and 
period  of  body  waves. 

The  influence  of  the  local  near  surface  geology  on  the 
amplitude  of  short  period  P  waves  was  observed  by  Fernandez 
(1963)  using  the  North  American  stations  of  the  Long  Range 
Seismic  Measurements  Program.  Similar  types  of  observa¬ 
tions  for  long-period  waves  by  Nuttll  (1964)  showed  that  the 
long  periods  are  not  so  sensitive  to  local  characteristics 
of  the  site.  The  apparent  angle  of  Incidence  and  its  changes 
with  frequency  were  the  object  of  observations  by  Nuttli 
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and  afhitmore  (i96l)and  by  Nuttli  (1964a). 

1.2  Theoretical  studies  of  tne  Influence  of  a  layered 
medium  on  the  spectrum  of  longitudinal  waves. 

The  series  of  observations  just  mentioned  were  accom¬ 
panied  by  theoretical  developments  of  the  subject. 

The  observations  of  Imamura  in  1929  attracted  the 
attention  of  several  Japanese  seismologists.  The  problem 
of  the  free  vibrations  of  the  layered  crust,  and  especially 
of  the  top  layer,  was  not  purely  academic.  In  a  country 
like  Japan  the  influence  of  the  ground  on  the  amplitudes  and 
accelerations  created  by  earthquakes  was  a  matter  of  public 
safety  and  of  practical  interest. 

Under  certain  simplifying  assumptions  the  theoretical 
problem  concerns  the  solution  of  the  wave  equation  which 
satisfies  the  boundary  conditions  of  a  system  of  horizontal 
layers.  The  layers  are  excited  by  an  incident  longitudinal 
wave  which  arrives  at  the  base  of  the  system  at  some  angle 
of  incidence.  The  boundary  conditions  for  the  solution  of 
the  problem  are  the  continuity  of  stress  and  displacement 
at  each  internal  interface  and  the  vanishing  of  the  stress 
at  the  free  surface.  The  numerical  calculations  required 
were  laborious  and  the  authors  simplified  the  problem  to  the 
simple  cases  of  one  or  two  layers  and  normal  angle  of  inci¬ 
dence  , 

This  is,  for  example,  the  approach  of  Sezawa  (I930)who  con¬ 
siders  a  one-layer  model  over  a  half-space  and  a  dllatationa.'. 


9 


pulse  propagated  vertically.  The  same  method  is  used  by 
Sezawa  and  Kanai  in  a  series  of  studies  (1932a,  b,  1937} 
including  two-layer  models,  obliquely  incident  waves  and 
harmonically  oscillating  sources.  The  theoretical  seismogram 
is  a  function  of  the  product  of  the  frequency  of  the  waves 
and  the  thickness  of  the  layers.  The  same  authors,  Sezawa 
and  Kanai,  in  another  series  of  articles  (1932c,  1934), 
studied  the  results  of  reflections  and  refractions  in  a 
stratified  medium. 

Kanai,  in  a  collection  of  papers  (1952,  1953a*  b),  cal¬ 
culates  the  resonance  curve  for  the  vertical  and  horizontal 
components.  His  results,  though  limited  and  partial  and  ob¬ 
tained  by  a  laborious  method,  are  exact  and  represent  a 
significant  contribution  to  the  transmission  problem  of 
seismic  waves  in  layered  media. 

Finally,  Kanai  and  Yoshizawa  (1956)  studied  the  same 
effect  for  three-layer  models  excited  not  only  by  infinite 
trains  of  harmonic  vibrations  but  also  by  finite  trains  of 
only  one  or  two,  or  at  most  a  few  complete  cycles.  Kanai, 
Tanaka  and  Yoshizawa  (1959)  extended  the  number  of  models 
and  made  similar  calculations  for  the  motion  at  the  free 
surface  and  at  the  Interfaces . 

In  the  meanwhile  a  new  approach  to  the  problem  was  de¬ 
veloped  by  Thomson  (1950),  and  by  Haskell  (1953).  The  meth¬ 
od  allows  the  calculation  of  the  solution  to  the  wave  equa¬ 
tion  for  any  number  of  layers  in  terms  of  a  four  by  four 
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matrix  whose  members  are  a  function  of  the  parameters  of 
each  layer,  the  angle  of  Incidence  and  wave  number  of  the 
wave.  The  method  Is  very  convenient  for  digital  computer 
calculations.  In  his  numerical  calculations  for  a  standard 
crustal  model  of  one-layer  over  a  semi-infinite  half-space 
Haskell  (1962)  gives  the  ratio  of  the  vertical  and  horizon¬ 
tal  components  of  ground  motion  to  the  components  of  incident 
sinusoidal  P-wave  at  the  bottom  of  the  layer,  as  a  function 
of  frequency  and  angle  of  incidence. 

The  work  of  Haskell  was  extended  to  different  layered 
models  by  Hannon  (1964a).  The  effects  of  the  layers  on  the 
vertical  and  on  the  horizontal  component  are  considered  as 
transfer  functions  and  the  influence  of  changes  in  the 
thicknesses  and  velocities  of  the  layers  is  carefully  exam¬ 
ined.  Theoretical  seismograms  have  been  obtained  by  Hannon 
(1964b)  by  a  Fourier  synthesis  of  the  transfer  functions 
computed  using  the  Haskell -Thoms on  matrices. 

Similar  calculations  for  different  crustal  models  were 
the  basis  for  the  crustal  studies  of  Phlnney  (1964)  who 
compared  theoretical  curves  for  the  ratio  of  the  vertical 
to  the  horizontal  transfer  functions  with  the  ratio  of  the 
spectra  of  vertical  and  horizontal  seismograms.  In  this  way 
Phinney  was  able  to  evaluate  the  crustal  structures  at  some 
localities.  Nuttli  (1964a)  used  the  transfer  function  to 
correct  the  apparent  angle  of  incidence  obtained  from  ooser- 
vations.  The  transfer  functions  for  the  S  wave  .iave  also 
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been  used  by  Nuttli  (1964b)  to  correct  the  angle  of  polar¬ 
ization  of  S  waves.  In  the  studies  of  focal  mechanisms 
from  the  spectrum  of  body  waves  Ben-Menahem,  Smith  and  Teng 
(1964)  correct  the  spectra  for  the  effects  of  the  crust  at 
the  station  site  using  the  transfer  functions. 

F  antly  Phinney  (1965),  to  facilitate  the  use  of  theo¬ 
retical  calculation,  integrated  numerically  the  solutions 
for  the  response  to  a  line  source  and  a  point  source  in 
elastic  media.  This  method  of  integration  gives  a  theoreti¬ 
cal  spectrum  of  the  first  seismic  arrival,  as  viewed  through 
an  exponentially  decaying  time  window. 

Russian  seismologists  have  shown  interest  in  the  prac¬ 
tical  applications  of  the  problem.  There  are  several  In¬ 
stances  in  which  they  have  used  the  influence  of  the  sedimen¬ 
tary  layers  on  the  body  waves  to  calculate  the  structure  of 
a  site  for  prospecting  purposes.  Savarensky  (1952)  studied 
the  changes  of  the  angle  of  emergency  with  frequency. 
Gamburtsev  (1954)  suggested  the  "seismic  frequency  sounding 
method"  to  Investigate  cross  sections  of  interest.  This 
method  consists  in  tue  comparison  of  observed  apparent  an¬ 
gles  of  Incidence  at  different  frequencies  with  theoretical 
values  calculated  using  layered  models.  Ivanova  (1959,1960) 
used  the  method  In  a  very  simple  way.  Observations  are  ob¬ 
tained  by  tuning  the  seismometer  at  different  frequencies 
and  observing  the  apparent  angle  of  Incidence  of  the  cor¬ 
responding  waves.  The  waves  in  question  are  generated  by 
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explosions  in  deep  wells.  The  results  are  compared  wit r 
theoretical  curves  obtained  using  a  method  developed  by 
Malinovskaya  (1959)* 

Halperln  (1962)  studied  the  changes  of  the  angle  of 
incidence  with  frequency  for  a  wave  passing  through  a  low 
velocity  layer.  An  interesting  theoretical  study  of  the 
effects  of  thin  layers  on  elastic  waves  is  given  by 
Pod'Yapol'skl  ( 1961 )  who  solves  the  Zoeppritz  equations 
for  complex  amplitudes  and  considers  the  effect  on  a  sinus¬ 
oidal  wave  multiply  reflected  in  a  thin  layer.  The  multiple 
reflection  introduces  the  frequency  effect  into  his  calcula¬ 
tions  * 

At  the  present  moment  it  can  be  said  that  our  theoreti¬ 
cal  knowledge  of  the  effect  of  a  layered  medium  on  longi¬ 
tudinal  waves  is  exact  and  complete.  The  Haskell  method 
offers  a  fast  technique  to  calculate  the  effect  of  the  lay¬ 
ers  on  the  vertical  and  horizontal  component  of  ground  mo¬ 
tion.  Nevertheless  the  theoretical  calculations  can  be  sys¬ 
tematized  in  a  better  way  in  order  to  group  together  into 
families  of  curves  models  with  identical  or  similar  response. 
For  this  purpose  use  may  be  made  of  a  dimensionless  parameter. 
Such  a  presentation  allows  a  better  application  of  the  theory 
to  observations  and  consequently  a  better  determination  of 
crustal  parameters.  Finally,  the  theory  should  be  applied 
to  more  refined  observations.  This  has  been  done  in  this 
study  using  earthquakes  of  intermediate  and  great  focal 
depth. 
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2.  Influence  of  the  Layer  Parameters  on  the 

Transfer  Functions 

The  Thoms on -Haskell  matrix  formulation  of  the  boundary 
problem  of  seismic  waves  transmitted  in  a  layered  medium 
offers  a  rapid  and  exact  calculation  of  the  corresponding 
transfer  functions.  Only  for  high  frequencies,  as  indi¬ 
cated  by  Dunkin  (1965)  is  it  possible  that  the  method  intro¬ 
duces  errors  due  to  numerical  calculations.  These  high 
frequencies  will  not  be  considered  here. 

A  detailed  exposition  of  the  method  has  been  given  by 
Haskell  (1962)  and  by  several  others  following  HasKell 
(e.g.,  Hannon,  1964)  and  will  not  be  repeated  here. 

In  order  to  investigate  crustal  structures  using  the 
data  of  crustal  transfer  functions  determined  from  the  ob¬ 
servation  of  P  waves  at  a  given  place.  It  would  be  desir¬ 
able  to  develop  an  inversion  program  whereby  from  the  ob¬ 
served  transfer  functions  one  might  determine  a  crustal 
model  which  best  fits,  in  a  least  squares  sense,  the  obser¬ 
vational  data.  To  this  end,  and  to  study  the  feasibility 
of  such  a  program,  it  is  necessary  to  understand  the  influ- 
enc  *hat  changes  of  the  parameters  of  the  layers  may  have 
on  the  appearance  of  the  transfer  function  curves.  In  order 
to  do  this  two  tests  have  been  applied.  In  the  first  a 
crustal  model  was  altered  by  a  small  fraction  either  in  the 
crustal  thickness  or  In  the  crustal  velocity  of  P  waves  and 
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the  results  compared  with  the  original  transfer  function. 

In  the  second  test  the  first  partial  derivatives  of  the 
transfer  functions  with  respect  to  the  thickness,  the  elas¬ 
tic  parameters  and  the  density  of  the  layers  were  calculated. 

2.1  Results  of  small  changes  of  the  parameters. 

The  transfer  functions  for  the  vertical  and  horizontal 
components  TW,  TU,  and  the  ratio  TW/TU  are  each  a  function 
of  several  parameters  here  Indicated  with  their  respective 
symbols  and  corresponding  to  the  ith  layer: 

h^  =  thickness  of  the  layer 

=  longitudinal  P  wave  velocity 
«  density 

-  shear  wave  velocity 

ln  =  the  angle  of  incidence  of  the  plane 

wave  at  the  base  of  the  layered  system. 

'  =  frequency  of  a  simple  harmonic  com¬ 

ponent  of  an  incident  wave. 

Each  combination  of  these  parameters  will  give  a  dif¬ 
ferent  value  of  the  transfer  functions.  For  a  given 
crustal  model  and  angle  of  incidence  the  values  of  the 
transfer  functions  versus  frequency  define  the  transfer 
function  curves . 

To  study  the  Influence  of  the  changes  of  the  above- 
mentioned  parameters  some  limitations  which  are  imposed  by 
nature  itself  may  simplify  the  problem. 


A  first  limitation  based  on  physical  properties  of  the 
rocks  of  the  crust  is  the  relationship  observed  between 
longitudina*  and  shear  waves.  It  is  commonly  accepted  that 
for  this  part  of  the  earth  the  ratio  ^  is  approximate¬ 

ly  constant  and  does  not  vary  greatly  from  V 3,  correspond¬ 
ing  to  a  value  of  0.23  for  the  Poisson's  ratio. 

Another  useful  relationship  based  on  observations  ex¬ 
ists  between  the  P  velocity  of  crustal  layers  and  the  den¬ 
sity.  According  to  Birch  (1964),  for  the  rocks  of  the  crust, 
the  density  increases  in  almost  a  linear  way  with  seismic 
velocities.  This  convenient  relation  will  be  used  in  our 
calculations. 

Transfer  functions,  calculated  using  the  FORTRAN  II 
program  listed  in  Appendix  I,  are  presented  in  the  graphs 
of  Figures  1  and  2.  In  Figure  1  the  curves  indicated  by 
the  asterisks,  *  ,  correspond  to  an  angle  of  incidence  at 
the  base  of  the  crust  of  35°  and  to  a  crustal  model  A  given 
by  the  parameters: 

P  S 

Thickness  Velocity  Velocity  Density 
h  (Km)  (Km/sec)  (Km/sec)  (gr/cc) 

Crustal  Layer  30  6.56  3.7872  2.64 

Mantle  00  8.2  4.7344  3.33 

The  curves  indicated  by  the  triangles,  £  ,  correspond  to 

model  B.  This  model  has  exactly  the  same  parameters  but 

with  the  thickness  of  the  crust  Increased  by  20#  to  36  Km. 

Examination  of  these  two  sets  of  curves  shows  that  the 
second  set  is  fore-shortened  or  "shrunk"  In  frequency 
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Modulus  and  phase  of  the  transfer-  functions  TW,  TU  and  TW/Tt)  for  the  crustal  models  A  and  C  and  for  rays  at  3C,  angle  of  Incidence. 
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dependence  but  Is  otherwise  exactl;,  the  same  in  amplitude 
and  shape  as  the  first.  This  is  obvious  not  only  for  the 
absolute  values  of  the  three  transfer  functions  TW,  TU,  and 
TW/TU,  but  also  for  the  phase  angle  presented  in  the  graph 
at  the  right  side  of  the  figure.  As  the  frequency  increases 
the  separation  of  similar  portions  or  the  curves  becomes 
greater.  This  Indicates  that  both  models  will  show  a  sim¬ 
ilar  spectrum  of  body  waves  at  low  frequencies  but  quite 
different  at  higher  frequencies.  The  shifting  of  the  fea¬ 
tures  of  the  curve  corresponding  to  the  altered  model  rela¬ 
tive  to  tho  original  one  is  exactly  20#  on  the  frequency 
scale.  If  the  plot  were  made  on  logarithmic  paper  this 
effect  would  be  shown  by  a  displacement  to  the  left  of  the 
second  curve;  otherwise  the  appearance  of  both  curves  would 
be  exactly  the  same. 

The  properties  of  the  matrices  used  to  calculate  the 
transfer  functions  already  suggested  this  factor,  as  was 
indicated  by  Haskell  (1953).  Similar  effects  are  observed 
when  multilayered  models  are  considered. 

In  a  similar  way  the  same  model  A  was  compared  with 
model  C.  This  model  is  exactly  as  model  A  but  the  P  veloc¬ 
ity  in  the  crust  was  incremented  by  2(X,  to  7.87  km/sec. 

The  results  of  this  change  are  presented  in  Figure  2.  The 
shear  velocity  and  the  density  are  also  changed  according 
to  the  relations  previously  indicated.  Examination  of  the 
results  shows  that  in  this  case  the  change  of  velocities 
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affects  not  only  the  frequencies  but  also  the  amplitudes  of 
the  curves. 

The  effect  on  the  frequencies  is  obvious,  especially 
in  the  phase  angle  curves .  In  this  case  the  frequency 
scale  is  expanded  by  20$. 

The  change  of  amplitude  of  the  transfer  function  os¬ 
cillations  was  expected,  since  the  contrast  of  velocities 
between  mantle  and  crust  is  smaller;  consequently  the  re¬ 
flection  coefficients  for  the  P  and  SV  reflected  within  the 
crust  are  smaller.  The  physical  interpretation  is  that  for 
small  contrast  of  velocities  the  greatest  part  of  secondary 
energy  is  transmitted  back  into  the  mantle. 

This  amplitude  effect  due  to  the  change  of  velocities 
influences  not  only  the  absolute  values  of  the  transfer 
functions  out  also  the  oscillations  of  the  phase  angle. 

The  same  two  tests  were  performed  for  the  same  models 
for  the  angle  of  incidence  of  10°.  The  results  are  similar 
and  are  presented  in  Figures  3  and  4. 

These  results  are  significant  since  they  show  that  the 
transfer  functions  have  a  special  relationship  to  frequency 
such  that  changes  of  the  parameters  displace  the  peaks  and 
troughs  of  the  curves  in  the  frequency  domain.  This  behavior 
of  the  transfer  functions  presents  difficulties  when  an  in¬ 
version  problem  is  attempted  in  terms  of  frequency. 

2.2  First  partial  derivatives. 

The  use  of  surface  wave  data  for  the  investigation  of 
earth  structure  has  been  improved  recently  by  the  develop- 
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merit  of  an  inversion  technique  based  on  a  least  squares 
curve  fitting  procedure  whereby  theoretical  group  and  phase 
velocity  dispersion  curves  obtained  from  selected  earth 
models  are  adjusted  so  as  to  best  fit  observed  dispersion 
curves.  Dorman  and  Ewing  (1962)  and  Brune  and  Dorman  (1963) 
applied  the  least  squares  technique  by  the  repeated  calcula¬ 
tion  of  partial  derivatives  of  phase  velocity  with  respect 
to  each  of  the  parameters  of  the  layers.  This  is  done  for 
each  iteration  of  the  calculations.  Anderson  (1964)  showed 
that  for  reasonable  variations  in  the  parameters  of  the 
layers,  the  partial  derivatives  of  the  phase  velocity  may 
be  considered  constant.  This  property  reduces  the  least 
squares  inversion  technique  calculations  and  has  been  used 
with  success,  for  example,  by  McEvilly  (1964)  to  evaluate 
the  crust-upper  mantle  structure  in  the  central  region  of 
the  United  States. 

In  order  to  use  the  body  wave  spectra  to  evaluate  the 
structure  below  the  point  of  observation  the  same  technique 
could  be  applied  to  obtain  the  best  fit  in  the  least  squares 
sense  between  observations  and  theoretical  curves.  As  in¬ 
dicated  above,  a  necessary  condition  to  use  the  simplified 
technique  of  Anderson  is  that  the  partial  derivatives  of 
the  transfer  functions  with  respect  to  the  parameters  of 
the  layers  remain  constant  for  reasonable  variations  of  the 
parameters.  This  property  was  investigated  by  the  direct 
calculation  of  the  partial  derivatives  of  the  transfer  func¬ 
tions  TW,  TU  and  the  ratio  TW/TU . 
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Using  the  notation  and  sign  conventions  of  Haskell 
(1953  1962)  and  Hannon  (1964a),,  the  transfer  functions 

for  the  horizontal  and  vertical  components  of  ground  mo¬ 
tions  at  a  given  frequency  are: 


r  2.  C  (vXl  -  j31 

p 

TWM  -  2c  GX.  -  Jj.) 

<y„  P 


(2-1) 


p  a,-  j«)  -  <x~ 


and  the  are  the  elements  of  the  matrix  product  obtained 

in  solving  the  boundary  value  problem  of  the  layer  system: 
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J  O-A-x 


-  CL 


(2-2 ) 
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k  nee  each  element  or  sub-matrix  am  of  (2-2)  depends 
on  the  parameters  of  only  one  layer  it  is  possible  to  fina 
the  first  partial  derivatives  of  J  and  D  with  respect  to  any 
of  the  independent  parameters  in  any  layer  m  by  finding  the 
partial  derivatives  of  the  element  matrix  am  and  substituting 
this  value  in  the  product  J. 

For  example,  the  partial  derivative  of  TU(uj)  with  re¬ 
spect  to  the  thickness  of  the  mth  layer  will  be: 


^TV(u))  =  i_c 
d  o/~. 


(2-3) 


and 


2lL 
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h-l 


*  *  * 


(2-4) 


In  order  to  be  able  to  differentiate  am  with  respect  to 
the  independent  parameters  of  the  layers,  each  element  of 
a^  should  be  expressed  in  terms  of  the  density  the  thick¬ 
ness  In*,  and  the  elastic  parameter  A*,  .  For  the  crust  it 
is  assumed  that  the  Lame  constants  have  identical  value, 
that  is  X  -  M 

Then,  for  example,  the  first  element  of  am  is: 
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(‘ 
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and  its  first  partial  derivative  with  respect  to  the  thick¬ 
ness  of  the  mth  layer  is: 


The  same  operation  may  be  performed  for  the  16  elements  of 
the  matrix  and  with  respect  to  the  independent  parameters 

.  These  calculations  have  been  combined 
in  a  FORTRAN  II  computer  program  listed  in  Appendix  (II), 
and  the  results  of  this  program  were  plotted  in  Figures  6  to 
8.  The  partial  derivatives  shown  in  these  figures  correspond 
to  the  3ame  crustal  model  A  studied  earlier  in  this  chapter 
by  variation  of  the  parameters. 

In  general  the  partial  derivatives  have  large  values 
even  for  low  frequencies.  The  values  of  the  derivative 
are  greater  for  the  ratio  of  the  vertical  and  horizontal 
components  TW/TU.  By  studying  the  curves  of  the  first  partial 
derivatives  it  may  be  noted  that  the  maxima  correspond  to 
the  maxima  of  the  corresponding  transfer  functions.  This 
could  be  expected  from  the  character  of  the  shifting  toward 
higher  or  lower  frequencies  of  the  peaks  of  the  transfer 
functions,  as  noted  above,  when  the  parameters  of  the  layers 
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first  partial  derivatives 
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Hodulu®  and  phase  of  the  flret  partial  derivatives  of  the  transfer  functions  TW,  TU  and  TW/TU  with  respect  to 
angle  of  Incidence. 
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are  changed.  This  shifting  will  affect  most  notably  the 
peaks  of  the  curves.  Also,  since  the  shift  increases  with 
frequency,  tne  partial  derivatives  increase  also  with  fre¬ 
quency. 

To  study  the  effect  of  changes  in  the  parameters  in  a 
multilayered  model  the  crustal  model  1  of  Hannon  (1964a) 
was  selected.  This  model  is  similar  to  the  crustal  struc¬ 
ture  of  the  central  United  States  and  its  parameters  arei 


Layer 

Thickness 

P 

Vel03ltv 

S 

Velocity 

Density 

Number 

{km) 

(km/sec) 

(km/sec) 

(gm/cm3) 

1 

1.0 

4.40 

2.50 

2.70 

2 

19.0 

6.20 

3.50 

2.00 

3 

18.0 

6.40 

3. TO 

2.90 

Half-space 

8.20 

4.00 

3.30 

The  modulus  and  the  phase  of  the  transfer  function  for  the 
vertical  and  horizontal  component  TW  and  TU  and  their  ratio 
TW/TU  are  given  in  Figure  9.  The  corresponding  first  par¬ 
tial  derivatives  with  respect  to  the  thickness.  Lame's  con¬ 
stant  ar.d  density  of  each  of  the  layers  are  presented  in 
Figure  10  to  18. 

These  calulations  show  that  thin  layers  have  smaller 
values  of  the  first  partial  derivative  than  do  thick  layers. 
This  is  explained,  again,  by  the  shifting  of  the  frequencies 
since  the  shifting  was  proportional  to  the  thickness  of  the 
layer  in  relation  to  the  total  thickness  of  the  system.  In 
the  case  of  a  thin  layer  this  shifting  is  small  and  as  such 
the  changes  of  the  transfer  functions  will  be  small. 
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The  periodical  character  of  the  first  partial  deriva¬ 
tives  may  be  explained  by  the  same  shifting  of  the  frequen¬ 
cies.  As  the  frequency  increases  the  peaks  are  displaced 
more  and  more  till  one  peak  is  displaced  to  the  frequency 
of  the  neighboring  peak  of  the  original  curve.  In  this 
event  the  first  partial  derivative  is  not  so  great.  This 
situation  corresponds  to  the  minima  of  the  first  partial 
derivative  curves. 

Among  the  different  parameters  the  changes  in  density 
affect  more  the  transfer  functions.  For  the  angle  of  inci¬ 
dence  of  20°  the  vertical  conqponent  transfer  function  is 
more  sensitive  to  changes  of  the  parameters  than  the  hori¬ 
zontal  component. 

To  test  the  stability  of  the  first  partial  derivatives 
when  small  changes  of  the  parameters  are  Introduced  In  the 
parameters  of  the  crust,  model  A  and  its  partial  derivatives 
were  compared  with  the  first  partial  derivatives  of  model  D. 
This  model  D  is  the  same  as  model  A  with  the  exception  of 
the  thickness  of  the  crust  which  was  Incremented  by  IQ %  to 
33  ton. 

The  numerical  results  showed  that  the  individual  values 
of  the  first  partial  derivatives  for  a  given  frequency  were 
quite  different  for  both  models.  Consequently  an  inversion 
process  based  on  the  assumption  that  the  first  partial 
derivatives  remain  constant  for  reasonable  changes  of  the 
layer  parameters  is  not  possible  if  the  transfer  functions 
are  calculated  as  a  function  of  frequency. 


3.  Analysis  of  Factors  Affecting  the  Character 
of  the  Transfer  Function  Curves 

Though  normal  mode  theory  and  the  Haskell -Thoms on 
matrix  development  is  the  more  convenient  for  numerical  cal¬ 
culations  of  the  transfer  functions  of  seismic  waves  in 
layered  media,  nevertheless  there  are  other  theoretical 
approaches  to  the  problem.  While  these  latter  are  more 
laborious  for  numerical  calculations,  in  some  respects  they 
give  clearer  Insight  into  the  physical  problem  and  into  the 
phenomena  taking  place  in  the  stratified  medium.  This 
better  understanding  of  the  problem  helps  to  evaluate  the 
influence  of  each  of  the  layer  parameters  in  the  transfer 
function  itself. 

In  this  section  we  present  two  different  approaches  to 
the  propagation  of  seismic  energy  in  layered  media  in  order 
to  understand  better  the  behavior  of  the  transfer  functions. 
The  first  one  is  based  on  considerations  of  constructive 
and  destructive  interference  of  the  primary  and  reflected 
waves  arriving  at  a  given  point  of  the  surface.  The  second 
approach  considers  the  Fourier  transform  of  a  delta  function 
pulse  multiply  reflected  and  refracted  in  the  layered  medium. 

3.1  Interference  patterns  of  multiply  reflected  seismic 
waves . 

The  effect  of  layered  media  on  an  Infinite  train  of 
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sinusoidal  longitudinal  waves  of  a  given  angular  frequency 
<jJ  nt'j  be  investigated  from  the  point  of  view  of  construc¬ 
tive  interference  of  the  direct  arrival  of  the  wave  at  the 
surface  end  the  simultaneous  arrivals  of  reflected  and  re- 
fractea  waves  at  the  same  point  of  the  surface.  The  result 
will  be  a  superposition  of  waves  of  the  same  frequency  but 
of  different  amplitude  and  phase. 

The  amplitude  of  the  reflected  and  refracted  rays  may 
be  found  by  solving  the  Zoeppritz  equations  (Zoeppritz, 

1919)  in  terms  of  the  angle  of  incidence  of  the  ray,  the 
type  of  conversion  taking  place,  and  the  contrast  of  veloc¬ 
ities  between  the  ti.o  media  at  the  interface.  For  a  ray 
suffering  several  reflections,  refractions,  and  conversions 
the  final  amplitude  will  be  the  product  of  the  appropriate 
coefficients  at  each  interface.  Graphical  and  numerical 
values  for  the  solutions  of  Zoeppritz’ s  equations  are  given 
by  several  authors  (e.g.,  Steinhart  and  Meyer,  1959;  McCamy, 
Meyer,  Smith,  1962;  Ccstain,  Cook  and  Algermissen.  1963). 

In  all  these  studies  the  frequency  of  the  wave  is 
neglected  since  the  authors  are  interested  only  in  the  am¬ 
plitude  or  transmission  coefficients  of  tne  interfaces. 

However,  the  resultant  amplitude  for  the  displace¬ 
ment  at  the  free  surface  of  a  layered  media  due  to  the 
superposition  of  vraves  of  the  same  frequency  but  of  vary¬ 
ing  amplitude  and  phase  is  a  function  of  both  the  ampli¬ 
tudes  and  the  phases  of  the  components. 
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3-11  The  single  Payer  case. 

Consider  a  simple  case.  Let  a  plane  wave  front  DE 
(Figure  19)  arrive  at  the  interface  of  a  single  layer 
model. 


Figure  19.  Seismic  transmission  in  a  single  layer  model. 

Let  R  be  a  point  at  the  surface  such  that  the  direct  P  and 
the  converted  SV  from  the  point  F  meet  together.  The  two 
rays,  PD  and  PF,  were  in  phase  at  the  wave  front  DE  but  they 
will  not  be  in  phase  when  they  arrive  at  R.  The  phase  dif¬ 
ference  is  a  function  of  the  travel  time  difference  At  of 
the  two  paths.  If  these  travel  times  are  TEFR  and  Tpp  their 


difference  will  be: 
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where  e^  is  the  angle  of  emergence  of  the  ray  below  the 
interface 

e1  and  f^  the  angles  of  emergence  in  the  layer 
for  the  P  and  SV  waves. 

Expression  (3-1)  can  be  simplified  using  the  relations 

-£L-  =  -JL  ;  (3-2) 

so  that 

At  *  (—  f,  -  <?] 

d,  U,  f  7 

Assuming  a  Poisson's  ratio  6“ *  0.25  as  an  approxi¬ 
mate  value  for  the  rocks  of  the  crust,  the  time  difference 
is  reduced  to 

At"  =  i^\/J  i-Ar^.  j!(  —  £(  ^ 

If  T  is  the  period  corresponding  to  the  angular  veloc¬ 
ity  b^the  phase  shift  £sy  at  R  between  the  F  and  SV  wave 
will  be  given  by: 
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where  f  is  the  natural  frequency  of  the  wave. 

In  a  similar  way  the  phase  shift  £^p  between  the 
direct  P  arrival  and  the  doubly  reflected  P  wave  (Figure 
10)  may  be  found. 


R 


Figure  20.  Triple  reflection  of  P  wave. 

In  this  case  the  phase  difference  is: 

^  _ku2  (3-4) 

Ok  « 

4,  -  •• ) 

For  the  simple  case  of  a  single  layer  the  only  possible 


combinations  of  rays  arriving  to  the  receiver  R,  are  P  and 
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SV  waves  that  were  reflected  one  or  more  times  at  the  surface 
and  at  the  interface . 

In  general  the  phase  difference  between  the  direct  P 
arrival  and  a  wave  that  travels  in  the  layer  n  times  as  P 
and  m  times  as  SV  is  given  by; 


l)  Phase  shift  due  to  n  times  as  P  =  Jl*  (n-ij  €, 


2)  Phase  shift  due  to  m  times  as  3V  =  |.  ^  l/j 

So  if  £np,  msv  is  the  phase  3hift  of  the  total  path: 


ef  +-  /w,Vl  hv.|, 


(3-5) 


Note  that  (3-3)  and  (3-4)  are  special  cases  of  the 


general  expression  (3-5). 


The  phase  difference  between  any  two  rays  arriving 
at  R  one  of  the  traveling  N  times  as  P  ard  N  times  as  SV 
and  the  other  n  times  as  P  and  m  times  as  SV  is  given  by; 


£  +•  2.  n  fi  =  ^|.  Jll  (N“  n  )  +■  l/J 

(*’-  ••• ) 


(3-6) 


If  a  continuous  train  of  sinusoidal  waves  arrives  at 
the  interface  the  amplitude  at  the  surface  will  be  controlled 
by  these  phase  differences  and  by  the  amplitudes  of  the  re¬ 
flections.  If  two  rays  are  In  phase  they  will  reinforce  each 
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other  to  give  large  amplitude;  if  they  are  180°  out  of  phase 
their  resultant  will  be  the  difference  of  the  amplitudes. 
When  all  the  possible  reflections  are  considered  and  the 
ground  motion  is  resolved  into  vertical  and  horizontal  com¬ 
ponents  and  the  effect  of  the  free  surface  is  considered, 
the  final  result  is  the  value  of  the  transfer  function  for 
the  vertical  and  for  the  horizontal  component  for  the  given 
frequency,  as  found  by  Haskell  (1962).  The  same  considera¬ 
tion  for  other  frequencies  will  give  the  transfer  function 
curves  in  terms  of  frequency.  For  the  single  layer  case 
the  problem  could  be  easily  solved  and  prepared  for  numeri¬ 
cal  calculations  since  the  amplitude  coefficients  of  the 
reflected  waves  become  too  small  after  a  few  reflections  and 
their  contribution  could  be  neglected.  The  problem  will  be 
more  complex  when  several  layers  are  considered.  This  is 
the  reason  why  for  numerical  calculations  it  is  more  con¬ 
venient  to  use  the  matrix  method  of  Haskell  and  Thomson. 
Nevertheless  the  expression  (3-6)  for  the  phase  shift  is 
significant  since  it  indicates  the  influence  on  the  transfer 
function  of  the  thickness  of  the  layer,  and  the  velocity  of 
the  P  wave,  and  the  angle  of  incidence. 

The  thickness  of  the  layer  h-^  is  a  common  factor  of 
all  terras  in  the  expressions  (3-5)  (3-6).  This  means  that 
the  same  phase  relation  could  be  obtained  for  a  model  with 
a  layer  of  different  thickness  provided  the  frequency  is 
changed  in  such  a  way  that  f.h^  remains  constant.  Since 
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the  thickness  of  the  layer  does  not  have  any  influence  on 
the  transnu  sion  coefficients  the  amplitudes  of  the  trans¬ 
fer  functions  will  be  the  same  and  the  values  will  simply 
be  displaced  in  the  frequency  domain.  This  was  numerically 
and  graphically  shown  in  the  graphs  of  Figure  1. 

The  influence  of  the  P  velocity  on  the  transfer  func¬ 
tions  is  different.  A  change  of  OC1#  affects  also  the  fre¬ 
quency  in  a  way  similar  to  the  thickness  though  in  the 
opposite  direction.  The  main  influence  of  the  change  in 
velocity,  however,  is  that  the  transmission  coefficients  de¬ 
pend  on  the  velocity  contrast  above  and  below  the  interface. 
The  curves,  accordingly,  will  not  only  be  displaced  in  the 
frequency  domain  but  the  amplitudes  will  also  be  affected 
by  a  constant  multiplier.  This  effect  may  be  examined  in 
Figure  2  where  the  transfer  functions  of  two  crustal  models 
were  compared. 

These  considerations  suggest  that  the  same  transfer 
functions  could  be  obtained  if  the  values  f^hj/Q^  and 
ol1^2  were  to  be  held  constant,  though  the  transfer  func¬ 
tions  so  obtained  would  not  correspond  to  the  same  values 
of  frequency.  If  only  2  changes,  the  amplitude  of 

the  transfer  function  will  change  but  not  the  periodical 
character  and  general  appearance  of  the  curve. 

This  analysis  of  the  transfer  functions  suggested  the 
calculations  in  terms  of  the  dimensionless  quantity  x  where: 

)f  =  f  (  sin  e*  4-  fa  Sin  f*) 


(3-7) 
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Since  ck^/ol  g  can  have  any  value  this  results  in  a 
family  of  curves  with  similar  periodicity  and  shape  but 
different  amplitude  of  the  oscillations.  A  family  of  curves 
is  obtained  for  each  angle  of  incidence  since  the  trans¬ 
mission  coefficients  change  with  angle  of  incidence.  The 
change  with  angle  of  Incidence  is  gradual.  This  makes  it 
possible  to  Interpolate  between  the  different  angles  of  in¬ 
cidence  considered  in  this  study.  The  set  of  these  curves 
is  presented  and  discussed  in  detail  in  the  next  section. 


3.12  The  multi-layered  case. 

The  same  considerations  that  were  used  for  the  one 
layer  problem  can  easily  be  extended  to  the  two  and  more 
layer  case. 

For  the  two-layer  case  Snell’s  law  shows  that 


Figure  21.  Multiple  reflection  and  conversions 

of  incident  P  waves. 
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Referring  to  Figure  21  from  D  to  R  the  P  wave  travels 
with  velocity^  Then  for  a  frequency  f,  the  phase  change 
£  pi  from  D  to  R  is 

u, 

(»=  o,<-) 

where  A  ^  ia  the  wave  length  of  P  In  the  first  layer.  If 
the  wave  travels  n^  times  as  P  in  this  layer,  the  total 
phase  change  in  this  portion  of  the  path  will  be 


£n  +  2>,o  =  -aa 


i 


p< 


oi 


2  rt  n '  -  P.  ni  h4 

u  <A ,  xiV\  e, 

(*>'-  o,  if ... ) 


Similarly,  for  a  wave  traveling  ng  times  as  P  in  the  second 
layer  the  corresponding  phase  change  £pm2,  will  be: 


v  l*z, 

et 

(*v-  4,... ) 

for  a  wave  traveling  m^  times  as  SV  in  layer  1,  the  phase 
is: 


pnL 


tn 


!■ 


n 


+  2rrn'  -  JnJiL 


0, 

(n'c  M,- • J 

and  for  an  SV,  mg  times  in  layer  2  the  phase  is: 


+  In 


"■  ■  j- 


*v%. 


^z 


(r\'o  o, -f,..  j 


(3-8) 
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The  phase  change  for  the  portion  of  the  path  EH  depends  on 
the  distance  GH,  that  is,  on  the  difference  of  the  projec¬ 
tions  onto  the  interface  of  the  crust  portions  of  the  two 
rays.  This  difference  for  two  rays,  one  traveling  as  P,  n^ 
times  in  layer  1  and  n2  times  in  layer  2,  and  as  SV,  and 
mg  times  respectively,  and  the  other  traveling  as  p,  N-^  and 


N2  times,  and  as  SV;  and  Mg  times,  is: 

n,)  c os  e4 

<-os  e* 

e, 

S  » 

4.  (hA,-  wi4)  cos  |< 

+.  K-^) 

COf  fl 

e, 

tv 

and  the  phase  £  EH 

corresponding  to  EH  is 

£*h  +  • 

l», 

£i  ,UN >■-*'-) 

i 

cos  ea 

0Lt 

'3-9) 

jlK 

-  m,J  c oS  Jl t 

4-  ^  —  *’v,» 

;  CM1/, 

(n'«  0,V...) 

The  total  phase  di^ferenae  £  for  the  two  rays,  from  (3-8) 
and  (3-9)  is: 


£  +  zWi  hi- 

Pm  e<  +  1/J  (A- 

V 

+  JlL 

cIt, 

“ 

(N»."  •%)**  *  \/J  (/Yla-rntJstv> 
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This  expression  for  the  phase  difference  of  all  the 
rays  arriving  at  the  point  R  suggests,  as  in  the  single 
layer  case,  that  it  may  be  advantageous  to  use  a  dimen¬ 
sionless  quantity  J  in  the  calculation  of  the  transfer 
functions  for  two-layer  models.  This  time  the  definition 
of  is: 


e.,  +  U3  sm 


et  +  i/I  f » w 


or  t-  lfz 


(3-10) 


if 


fr'l  £'  +  V/ 3  il'v-k 


When  calculations  are  performed  with  models  of  exact¬ 
ly  the  same  values  for  the  h-^/n^  and^/Wg  ratios,  the 
transfer  functions  have  similar  form  provided  the  angles  of 
emergence  are  the  same  through  the  layers  of  the  crust. 

The  amplitude  of  the  oscillations  depends  on  the  velocity 
contrast  of  the  two  layers  with  the  mantle.  If  this  con¬ 
trast  is  kept  constant  the  same  numerical  values,  for  a 
given  t  »  are  obtained  no  matter  what  the  individual  values 
of  the  thicknesses  and  velocities}. 

The  same  considerations  apply  in  obtaining  a  dimen¬ 
sionless  parameter,  ^  ,  for  the  case  of  an  n-layered  model. 


In  this  case 


Consequently,  if  all  layer  thicknesses  are  multiplied  by  a 
constant,  R,  the  same  results  will  be  obtained  for  equal 
values  of  f  .  In  terms  of  frequency  the  change  will  be  in¬ 
versely  proportional  to  R.  If  all  the  velocities  including 
the  velocity  of  the  semi-infinite  medium  are  multiplied  by 
a  constant,  the  same  transfer  function  is  obtained  for  the 
values  of  if  since  the  angles  of  emergence,  e±  and  f1# 
through  the  layers  remain  constant.  If  the  velocity  con¬ 
trast  with  the  mantle  is  changed  while  the  velocity  con¬ 
trasts  within  the  layers  remain  constant,  families  of  .  irves 
are  obtained  with  similar  form  but  different  amplitude.  The 
more  layers  are  present  the  more  combinations  are  possible 
and  the  number  of  families  Increases  greatly. 

3.2  The  Fourier  transform  of  a  pulse  and  its  reflections. 

The  problem  of  the  transfer  function  of  layered  media 
for  seismic  waves  could  be  considered  as  the  ratio  of  the 
Fourier  transforms  of  the  Incident  pulse  at  the  bottom  of 
the  system  and  the  Fourier  transform  of  the  record  obtained 
at  the  surface  of  the  layers. 
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For  simplicity  sake  it  is  convenient  to  choose  a  unit 
delta  function  for  the  incident  pulse  s^nce  the  Fourier 
transform  of  this  function  has  unity  as  modulus  and  zero 
phase  for  all  the  frequencies  (see  Figure  22a, b). 


(h 


Figure  22.  a)  Modulus  of  the  delta  function 

Fourier  transform 
b)  Phase 


At  the  surface  the  record  will  consist  of  a  first  ar¬ 
rival  followed  by  several  reflections  and  conversions  of 
energy  arriving  with  certain  time  delays  with  respect  to 
the  initial  P.  Let  us  assume  that  no  critical  values  of 
the  angle  of  incidence  are  reached  in  any  of  the  layers; 
consequently  no  total  reflection  takes  place,  with  the  ap¬ 
propriate  phase  shift,  depending  on  the  angle  of  Incidence. 
Then  the  only  possible  phase  shifts  are  0  or  J\  ,  which  on 
the  record  will  be  shown  as  pulses  of  the  same  sign  as  that 
of  the  direct  or  initial  P  or  pulses  of  the  opposite  sign. 
The  vertical  component  record,  neglecting  smali  absorptions 
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that  could  take  place  in  the  crust,  will  look  like  Figure 
and  its  time  domain  function  will  be  of  the  type: 

J(t)  «  CL  f(k)  +  k  f  C  +...  (3 

where  the  /T>  are  the  time  lags,  and  a,  b,  ...  are  Zoep- 
pritz  coefficients. 
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Figure  23.  Typical  record  of  pulse  multiply  reflected. 

Considering  only  the  first  arrival  and  the  first  sec¬ 
ondary  wave  let  us  find  the  corresponding  Fourier  transform 
in  the  frequency  domain.  Let 

|,(fc)  =  a.  Set)  t  \>cf[t-%J 


Then 


F(wJ  -  a-+-p=  /  Sit-K)  £  d  t~ 

'  ten  J 


-JO 
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^et  us  change  the  variable: 

t'=  t  - 


Then: 


R  (»)  • 


dfc' 


j  _ i W  j 

'  (X4-b  e  ^  O.+.  b 


This  is  a  complex  quantity.  Its  modulus  is: 

J  If  ,  ~k  .2.  2.  ) 

PJ  (w>j|=i  'Jcc  4-  b  cos(tu1^J  ■+-  b 

-  £  +  b  coiX(co^J  4-  ■+' 


ip  Sin  (u/>jJ 


->z 


a.“ +•  b  +  2-cc  b  co^s 


.(3-13) 

•1 


and  the  phase  is 


-i 


y>  =  1-a n  _  b  sm 

dX.  4-  b  cqs  (yj%) 


(3-14) 


The  transfer  functions  presented  by  Hannon  (1964a)  are  the 
ratio  of  the  vertical  and  horizontal  component  at  the  sur¬ 
face,  to  the  total  amplitude  at  the  bottom  of  the  layers. 
Our  expression  for  F-^co)  represents  this  same  ratio  of  the 
surface  components  to  the  total  amplitude  at  the  bottom 


WliihllftllM 
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since  the  spectrum  of  the  amplitude  for  the  delta  function 
has  unity  as  modulus  and  zero  as  phase. 

For  the  simple  case  just  considered  of  a  single  re¬ 
flection  the  modulus  of  the  transfer  function  consists  of 
the  square  root  of  a  function  that  oscillates  sinusoidally 
about  the  line  a2  +  b2,  with  amplitude  2ab  and  with  period 
determinable  from  the  quantity  .  For  large  values  of 

the  coefficient  b  the  amplitude  of  the  oscillation  will  be 
large.  This  would  correspond  to  the  converted  SV  wave  or 
to  the  first  reflected  P  wave  since  the  amplitude  of  the 
reflections  and  refractions  decays  rapidly  after  a  few 
partitions  of  energy  at  interfaces  nr-  the  free  surface.  In 
general,  the  larger  the  contrast  between  the  interfaces  the 
larger  is  the  coefficient  for  the  reflections,  and  the  lar¬ 
ger  the. oscillations  of  the  transfer  functions.  Further¬ 
more,  the  frequency  of  the  oscillations  will  be  more  rapid 
with  larger  ^  .  This  corresponds  to  long  time  lags  be¬ 
tween  first  arrival  and  reflections,  and  shows  that  in  gen¬ 
eral  very  fast  oscillations  of  the  transfer  function  will 
have  small  amplitude  and  long-period  components  will  have 
larger  amplitude. 

The  same  argument  is  easily  extended  to  the  considera¬ 
tion  of  several  reflections.  For  a  record  of  the  type  of 
Figure  23  on  which  multiple  reflections  of  the  primary 
pulse  are  shown,  the  time  function  of  the  record  will  be 
as  in  equation  (3-12)  and  its  Fourier  transform  in  the 
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frequency  domain  has  as  modulus 
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(3-15) 


This  modulus  is  the  square  root  of  a  function  composed 
of  the  constant  value  (a2  4  b2  4  c2  4  . . . . )  plus  a  series 
of  sinusoidal  components  of  different  period  and  amplitude. 
In  practical  calculations  only  a  few  of  these  sinusoidal 
components  should  be  considered  since  they  become  very  small 
after  a  few  reflections  and  refractions. 

The  phase  angle  for  the  same  transfer  function  is: 


^  _  ier^'  bfiw  4  c  4- . .  . 

CL  4  b  coj  (lu^+  c  (3-16) 

This  phase  angle  is  zero  for  zero  frequency  and  indicates 
that  long  wave  lengths  are  not  affected  by  the  layers.  As 
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the  frequency  Increases  the  numerator  of  the  expression 
increases  and  the  denominator  decreases,  ao  the  phase 
angle  increases.  But  again  the  value  of  the  phase  angle 
is  controlled  by  the  sinusoidal  components  with  the  larger 
amplitude  terms.  This  can  be  observed  in  the  graphs  of 
the  transfer  functions  presented  in  the  next  chapter. 

Since  the  periodicity  of  the  modulus  and  phase  of  the 
transfer  functions  depends  on  the  time  lags  between  first 
arrival  and  later  reflections  and  refractions,  it  is  con¬ 
venient  to  calculate  these  time  lags  as  a  function  of  the 
layer  parameters. 

3.21  Time  lags  between  primary  arrival  and  multiple 
reflections , 

The  first  arrival  o*  energy  after  the  direct  P  will  be, 
in  the  case  of  the  one-layer  model,  the  converted  SV.  From 
Figure  19  and  equation  (3-3)  this  time  interval  A  t  is 
(assuming  0.25 ) 

rsv  =  k-fi/S"  si*  l  -  si*  ej 

<  V  0  (3-17) 

The  time  interval  fo.»  the  wave  reflected  as  P  at  the  free 
surface  and  as  P  again  at  the  bottom  of  the  layer  (see 
Figure  20  and  equation  (3-4))  is  given  by 

rJP  =  k.  2  sin  e, 


(3-18) 
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In  general  the  time  interval  between  the  first  direct 
P  and  a  ray  reflected  n  times  as  P  and  m  times  as  SY,  where 
n+m  *  2r,  r  =  C,  1,  2,  ..  is: 

(3-19) 

and  the  period  of  the  transfer  function  component  will  be, 
from  (3-1 p),  a  function  of 

uu 

Tills  expression  indicates,  again,  that  if  for  different  one- 
layer  crustal  models  the  quantities  sin  e-,  and  sin  fn  remain 
constant  the  same  periodicity  or  character  y  .  the  transfer 
function  will  be  obtained  when  the  results  are  calculated 
in  terms  of  the  dimensionless  parameter 

-  J-  [  s  *  n  (*4  +  V3  s  i  n  J 

The  amplitude  of  the  oscillations  will  depend  on  the  veloc¬ 
ity  contrast  between  the  layer  and  the  semi-infinite  half¬ 
space. 

3.22  The  multi-layered  case. 

The  problem  could  be  extended  to  consider  multi-layered 
models.  For  the  two-layer  case  let  us  consider  the  time  lag 
between  the  first  direct  P  arrival  and  a  ray  reflected,  re¬ 
spectively: 
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and  m  times  as  P  and  as  SV  in  the  first  layer, 
n 2  and  mg  times  as  P  and  as  SV  in  the  second  layer, 
then  the  time  lag  of  the  reflected  wave  is: 


/V 


1  *\P,  = 


U. 


*  i 


(h,-  i)  Sin  <?,+  \/j  ^ 
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l)  I  I  W  V5  Witxm  f3 


(3-20) 


and  the  sinusoidal  component  of  the  transfer  function  in¬ 
troduced  by  tills  reflection  for  a  constant  time  lag  will 
oscillate  in  the  uJ  dimension  with  a  period 


(3-21) 


Since  the  transfer  functions  are  plotted  against  the  param¬ 
eter  o  it  is  convenient  to  express  this  period  in  terms 
cf  ^ units.  Then  using  (3-11),  (3-13)  and  (3-21)  the 
period  of  the  sinusoidal  component  will  be 

nr)-  T (-)i  4  ii 


(3-22) 


For  example,  for  the  first  SV  arrival  in  a  single  layer 
crust  the  time  lag  Tsv  given  by  (3-17)  was: 

T,„  =  J£l  ( #3  sm  -  si  we*") 
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The  period  of  the  transfer  function  component  introduced 
by  the  first  SV  reflection  measured  in  ^  units  will  be: 

T(y)  =  (sin  e,  t  l/J  Jin  |l)  (3-23) 

V 3  -  Jin 

Similarly  the  period  of  the  transfer  function  components 
for  the  case  of  multilayered  models  is  given  by: 

TU)  r  —  v/«"  ^ 

in  _  (3-24) 

r 

where  0"  is  given  by  (3-20)  in  terms  of  the  thicknesses  of 
the  layers,  h^,  and  the  P  velocities  in  these  layers,  ctL  . 

This  analysis  of  the  sinusoidal  components  of  the 
transfer  function  is  used  to  classify  the  multilayered 
models  in  groups  with  identical  or  similar  transfer  func¬ 
tion  curves. 

Since  the  thickness  of  the  layers  does  not  affect  the 
amplitude  of  the  components  but  only  their  periodicity  ac¬ 
cording  to  (3-24),  identical  transfer  curves  will  be  ob¬ 
tained  when  all  the  thicknesses  of  the  layers  are  multi¬ 
plied  by  the  same  factor.  Thus  the  significant  effect  of 
the  thicknesses  of  the  layers  of  the  models  is  not  their 
absolute  value  but  their  ratio. 

Charge  in  the  velocities  of  longitudinal  waves  in  the 
layers  may  change  the  amplitude  of  the  components  of  the 
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transfer  functions.  Nevertheless  if  their  ratio  is  kept  con¬ 
stant,  the  amplitudes  will  remain  constant  and  also  the 
periodicity  of  the  components,  according  to  (3-24).  Thus 
the  significant  effect  of  the  velocities  is  also  their  ratio. 
Nevertheless  it  is  possible  to  keep  constant  the  periodicity 
of  the  components  and  at  the  same  time  to  change  their 
anplitude.  Since  in  (3-24)  the  period  does  not  depend  on 
the  velocity  of  longitudinal  waves  in  the  mantle  it  is  pos¬ 
sible  so  change  this  velocity  without  affecting  the  periodic¬ 
ity.  But  changing  the  contrasts  of  velocity  between  tho 
mantle  and  the  crust  layers  will  affect  the  amplitude  of 
the  sinusoidal  components  of  the  transfer  function  (3-15) 
(3-16). 

These  observations  are  the  basis  for  the  evaluation  of 
the  transfer  functions  lor  two  and  more  layer  models.  By 
changing  the  contrast  of  velocities  between  the  mantle  and 
the  layers  but  keeping  constant  the  relationship  between 
the  velocities  of  the  layers  of  the  crust,  families  of 
curves  of  identical  shape  but  of  different  amplitude  are 
obtained.  As  before,  these  allow  interpolation  for  inter¬ 
mediate  models. 
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4.  Master  Curves  for  the  Transmission  oi‘  Seismic 
Longitudinal  Waves  in  Layered  Media 

The  theore  leal  considerations  of  the  preceding  section 
prepare  the  way  for  a  presentation  of  the  transfer  functions 
for  different  crustal  models  in  a  methodical  way. 

In  the  present  section  we  shall  first  present  and  dis¬ 
cuss  the  complete  set  of  single-layer  crustal  models.  Then 
we  will  discuss  the  case  of  two-layer  models,  and  several 
families  of  this  type  will  be  presented. 

Unless  otherwise  indicated  the  following  general  assump 
tions  will  be  presumed  to  hold: 

1)  The  layers  are  horizontally  stratified. 

2)  Poisson* s  ratio  in  the  cru3t  is  equal  to  0.25. 

3) -  The  density  changes  linearly  with  the  P  velocity. 

4)  The  apparent  surface  velocity  must  be  greater 
than  the  P  velocity  in  the  top  layer  and  greater 
than  the  S  velocity  in  any  layer. 

The  second  assumption  is  based  on  numerous  measurements 
made  on  the  rocks  of  the  crust  and  upper  mantle.  The  third 
assumption,  as  noted  in  section  3#  is  based  on  experimental 
studies  of  the  changes  of  density  in  the  crust  with  seit- 
mic  velocities  by  Birch  (1964)  and  Nafe  and  Drake  as  cl ltd 
by  Talwani  et  al.  (1959).  bmull  departures  from  these  two 
assumptions  such  as  may  present  in  nature  do  not  altei 
the  values  of  the  transfer  function  in  an  appreciable  wa  -• . 
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Assumption  4  is  placed  in  order  to  avoid  the  total  re¬ 
flection  which  will  take  place,  even  for  the  SV  converted 
ray,  if  the  assumption  is  violated. 

In  order  to  study  the  transfer  functions  the  FORTRAN  II 
Computer  Program  of  Hannon  (1964a)  was  modified  to  calcu¬ 
late  these  functions  in  terms  of  the  dimensionless  param¬ 
eter  t  as  defined  in  section  3  (see  Appendix  I).  The  pro¬ 
gram  can  also  be  used  to  calculate  the  functions  in  terms  of 
frequency  with  the  use  of  an  optional  control.  The  results 
are  plotted  on  logarithmic  scale. 

In  discussing  the  transfer  function  curves  we  shall 
make  use  of  the  terms,  periodicity  and  amplitude.  Perio¬ 
dicity  refers  to  the  oscillatory  character  of  the  curves 
with  their  peaks  and  troughs.  The  term  amplitude  refers 
to  the  maximum  values  of  the  transfer  functions  over  these 
oscillations. 

4.1  One-layer  models. 

The  simplest  crustal  model  is  a  single  layer  over  a 
half-space.  In  this  layer  the  velocity  for  the  P  and  S 
waves  may  be  considered  to  correspond  to  the  average  veloc¬ 
ity  of  the  P  and  S  waves  in  the  crust  of  the  real  earth, 
whether  the  crust  may  actually  consist  of  several  layers  or 
whether  there  may  exist  velocity  gradients  in  the  crust. 

The  thickness  of  the  layer  is  equal  to  the  total  thickness 
of  the  real  crust.  Though  this  is  a  simple  model,  never¬ 
theless  fr^om  the  theory  of  transfer  functions  given  in  the 
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last  chapter  It  may  be  expected  that  the  small  velocity 
contrasts  which  may  be  present  within  the  crust  would  in¬ 
troduce  only  small  amplitude  components  into  the  transfer 
function  which  will  not  al  ,er  substantially  the  general 
appearance  of  the  transfer  function.  The  largest  contrast 
of  velocity  which  might  be  expected  within  the  crust  would 
be  found  near  the  top  cf  the  crust  between  the  sedimentary 
layers  and  the  basement  rock,  but  these  layers  are  thin  in 
comparison  to  the  total  thickness  of  the  crust  and  conse¬ 
quently  their  influence,  though  important  in  terms  of  am- 
pj'tude,  will  be  spread  over  a  long-period  component  which 
not  affect  to  any  great  extent  the  appearance  of  the 
low  frequency  region.  These  theoretical  suggestions  are 
confirmed  by  actual  calculation  wnen  the  transfer  functions 
of  the  two-layered  models  are  compared  with  the  correspond¬ 
ing  ones  of  the  one-layer  models .  The  main  features  of  both 
are  common  and  correspond  to  the  large  velocity  contrast 
existing  between  the  mantle  and  the  crust. 

The  set  of  one  layer  models  (Figures  24  to  36)  is  com¬ 
plete  in  the  sense  that  for  each  angle  of  incidence  a  single 
family  of  curves  exists  with  identical  periodicity  which 
represents  all  possible  crustal  thicknesses  and  all  possible 
crust-mantlo  velocity  contrasts.  For  example,  in  Figure  30, 
corresponding  to  an  angle  cf  incidence  of  30°  at  the  bottom 
of  the  crust,  all  possible  models  of  different  thickness 
will  give  the  same  type  of  curve;  only  the  amplitude  will 
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Figure Tu  :  Modulus  and  phase  of  the  transfer  functions  TW,  TU  snd  TW/TU  for  one  layer  crustal  model  and  for  1  an*rle  of  incidence. 
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Figure  26 :  Modulus  ar.d  phase  of  the  transfer  functions  TW,  TU  and  TW/YU  for  one  layer  crustal  models  and  for  10  angle  of  incldene 
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Figure  36  !  Modulus  and  pheae  of  the  transfer  functions  TV,  TU  and  TV/TU  for  one  layer  crustal  mode’a  and  for  60°  angle  of  ; 


change  as  the  ratio  between  the  velocities  of  the  crust  and 
mantle  changes.  Exact  curves  are  given  for  the  ratios 

eQual  0.9  and  0.7  and  Intermediate  values  may  be 
Interpolated . 

In  the  non-realistic  case  in  which  the  P  velocity  in 

the  crust  is  greater  than  in  the  mantle,  the  sign  of  the 

amplitude  is  reversed.  For  illustrative  purposes  one  such 

o 

case  was  calculated  for  the  angle  of  incidence  of  25  and 
for  the  ratio  &^/o(  2  =  1.1  and  is  presented  in  Figure  2y. 

In  each  figure  the  modulus  and  phase  of  the  vertical 
and  horizontal  component  transfer  functions  are  given. 

These  expressions  are  theoretically  important,  but  in  practi 
the  spectrum  of  the  apparent  angle  of  emergence  is  the  one 
which  can  be  obtained  most  directly  from  observations  and 
special  attention  is  directed  to  this  portion  of  each  figure 
The  vertical  and  horizontal  component  are  given  in  order  to 
understand  better  the  features  of  the  apparent  angle  of 
emergence . 

In  these  figures  the  values  of  the  vertical  and  hori¬ 
zontal  component  are  normalized  so  as  to  represent  the  ratio 
of  the  surface  component  to  the  corresponding  component  of 
the  incident  motion  at  the  base  of  the  layered  system.  In 
this  way  for  infinitely  long  wave  lengths  the  value  of  the 
ratio  is  always  unity.  Independent  of  the  angle  of  incidence 
of  the  ray.  The  transfer  function  for  the  apparent  angle 
of  emergence  is  similarly  normalized  to  the  v.  ue  of  the 
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tangent  of  the  apparent  angle  of  emergence  corresponding  to 
Infinitely  long  wave  lengths;  this  corresponds  to  the  tan¬ 
gent  of  the  apparent  angle  of  emergence  at  the  bottom  of 
the  layers,  which  may  be  deduced  from  simple  ray  theory 
without  frequency  considerations  and  given  by: 

2.  (.0$ 1  t  =.  3  (<  -  siv\  e.) 

where  e  is  the  real  angle  cf  emergence  and  e  the  apparent 
one. 

4.11  The  vertical  component  transfer  function. 

Since  the  transfer  function  depends  on  several  varia¬ 
bles  such  as  the  thickness  of  the  layer,  velocities  of  the 
seismic  waves,  angle  of  incidence  of  the  ray,  and  the  den¬ 
sities,  of  the  materials,  it  Is  convenient  to  remark  the 
effect  of  these  parameters  and  the  effect  of  their  changes 
on  the  transfer  functions  The  remarks  which  follow  may  be 
verified  by  careful  examination  of  the  figures. 

l)  For  small  angles  of  incidence  such  as  5°  (Figure  25) 
the  main  contribution  to  the  vertical  component 
after  the  initial  P  arrival  pertains  to  the  P  wave 
reflected  once  at  the  surface  and  again  at  the 
crust-mantle  boundary.  The  amplitude  of  this  re 
flection  depends  on  the  contrast  of  velocity  be¬ 
tween  crust  and  mantle  and  in  each  case  could  be 
obtained  from  the  graphical  solutions  of  Zo?ppritz's 
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equation  mentioned  in  section  3.  In  our  graphs 
this  contribution  is  easily  detected  and  corresponds 
to  the  amplitude  of  the  sinusoidal  fluctuation  of 
the  vertical  component  transfer  function.  Since 
for  small  angles  of  incidence  the  reflected  P  is 
almost  the  only  reflection  of  importance,  it  occurs 
without  interference  from  other  reflected  phases. 
This  is  also  the  reason  why  the  amplitudes  of  all 
the  peaks  are  the  same. 

2)  For  small  angles  of  incidence  the  periodicity  of 
the  contribution  of  the  P  reflection  is  easily  cal¬ 
culated  from  the  time  lag  with  respect  to  the 
direct  P  arrival.  The  "period"  of  the  oscillation 
may  be  obtained  from  formula  (3-22).  For  an  angle 
of  incidence  of  5°  this  "period"  or  peak-spacing 


is: 


T  (y\  -  ytrie*  +  fL 
*"  XsmCj 


1.36, 


measured  in  y  units .  This  is  the  period  of  the 
predominant  large  sinusoidal  oscillation  of  the 
transfer  function  observable  in  Figure  23. 

3)  The  larger  the  contrast  of  velocities  the  larger 
is  the  amplitude  of  the  oscillations.  This  effect 
obeys  an  almost  linear  relation  which  allows  the 
interpolation  for  intermediate  values  of  the  con¬ 
trast  of  velocities.  If  the  ratio  of  velocities  js 
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reversed  so  that  the  larger  values  correspond  to  the 
crust  (see  Figure  29),  the  transfer  function  obtained 
is  symmetric  to  the  normal  case,  the  axis  of  symme¬ 
try  being  the  line  TW  =  1.  This  line  corresponds  to 
the  transfer  function  of  the  model  with  a  crust  of 
equal  P  velocity  to  that  of  the  mantle.  This  symme¬ 
try  of  the  transfer  functions  for  reversals  of  veloc¬ 
ity  implies  that  we  may  expect  the  transfer  functions 
to  be  sensitive  to  low  velocity  channels  at  depth. 

4)  The  influence  of  the  converted  SV  on  the  vertical 
component  is  very  small  at  small  angles  of  incidence. 
This  influence  increases  as  the  angle  of  incidence 
increases  and  is  quite  obvious  at  20°  or  25°  where 
small  oscillations  of  the  transfer  function  are 
superimposed  on  the  large  ones.  Nevertheless  the 
calculation  of  the  corresponding  "period"  shows  that 
these  short-period  oscillations  correspond  not  to 
the  initial  SV  conversion  but  to  combinations  of  P 
and  SV  reflections.  For  an  angle  of  incidence  of 
25°  the  different  "periods"  corresponding  to  these 
waves  deduced  from  their  time  lags  and  equation 
(3-22)  are: 

Ts0Tj  *3.3  ^units 

t1P,2S(*)=  -77  f- units 

1.0  units 

Thus  in  Figure  29  short  period  oscillations  of  the 
vertical  transfer  function  with  peaks  at  1.3  and 
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1.8  correspond  rather  to  the  multireflected  waves 
than  to  the  long-period  contribution  of  the  direct 
SV  converted  wave. 

5)  As  the  angle  of  incidence  increases  not  only  does 
the  vertical  component  of  the  SV  wave  increase  rela¬ 
tive  to  the  horizontal  component,  but  also  the  con¬ 
version  factor  of  P  into  SV  increases.  This  ex¬ 
plains  the  influence  of  SV  conversions  as  the  angle 
of  incidence  increases. 

4.12  The  horizontal  component  transfer  function. 

Due  to  the  influence  of  the  SV  conversions  from  the 
original  incident  P,  the  transfer  function  of  the  horizontal 
component  is  mere  complex  than  that  of  the  vertical  compo¬ 
nent. 

The  behavior  of  the  horizontal  component  has  a  special 
influence  on  the  value  of  the  tangent  of  the  apparent  angle 
of  emergence  iice  the  peaks  of  this  curve  correspond  to 
the  troughs  of  the  horizontal  transfer  function. 

Some  of  the  main  characteristics  of  this  component 
transfer  function  curves  are: 

l)  The  main  contribution  to  the  character  of  the 
curve  corresponds  to  the  SV  converted  wave,  at 
least  for  small  angles  of  incidence.  This  con¬ 
version,  nevertheless,  has  small  transmission 
coefficients,  as  indicated  by  the  graphical 
display  of  the  computations  of  MeCamy,  et  al.. 


(1962)  mentioned  previously. 

The  period  corresponding  to  the  first  SV  con¬ 
version  is  given  by  the  formula  (3-22) 

(sin  e.  +  i/3  sin  f, ) 

T  =!  1  1 

s  v  . '  rrr - 

(  V  3  sin  fj  -  sin  e^ ) 

For  an  angle  of  Incidence  of  25°  the  period  is 
3.3  f't  and  this  is  the  period  of  the  large 
oscillation  evident  on  the  graph  of  the  trans¬ 
fer  function  for  this  angle  of  incidence  (see 
Figure  29).  The  periodicity  of  the  transfer 
function  of  the  horizontal  component  is  not  so 
clear  as  was  that  of  the  vertical  component  since 
there  are  present  many  oscillations  of  almost 
the  same  amplitude. 

The  short-period  oscillations  of  the  transfer 

function  correspond  to  the  horizontal  components 

of  the  reflections  of  the  P-SV  type.  For  example, 

the  corresponding  periods  of  these  reflections 

o 

for  the  case  of  an  angle  of  incidence  of  25  are: 

T  1  gj*)-  .77  f units 

Tp,aCr)~  1'°  ^units 

These  short-period  components  are  very  clear  in 
the  graphs  (Figure  29). 

As  was  the  case  for  the  vertical  component,  the 
transfer  functions  of  larger  amplitude  correspond 
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to  the  cases  of  the  larger  contrasts  of  crust 
mantle  f  velocity.  This  is  obvious  not  only  in 
the  curves  of  the  modulus  but  also  in  the  phase 
angle  graphs. 

4 . 13  Curves  for  the  tangent  of  the  apparent  angle  of 
emergence . 

For  practical  purposes  the  curve  for  the  apparent 
angle  of  emergence  is  the  more  important  since  it  can  be 
obtained  directly  from  the  seismograms. 

The  main  features  of  the  curve  depend  on  the  angle  of 
incidence  of  the  P  wave.  This  angle  is  a  function  of  the 
epicentral  distance  and  of  the  depth  of  the  focus. 

Some  characteristics  of  these  curves  and  their  de¬ 
pendence  on  the  parameters  of  the  crust  may  be  indicated: 

1)  For  small  angles  of  incidence  of  the  ray  corres¬ 
ponding  to  large  epicentral  distances,  the 
spectrum  of  the  apparent  angle  of  emergence 
oscillates  rapidly  between  values  above  and 
below  unity.  Unity  corresponds  to  the  tangent 
cf  the  apparent  angle  of  emergence  calculated 
from  simple  ray  theory. 

At  shorter  epicentral  distances  the  apparent 
angle  of  emergence  is  more  stable  and  in  general 
remains  above  unity. 

2)  The  preceding  remarks  apply  only  to  a  normal  type 
of  crust  where  the  P  velocity  of  the  crust  is 


lower  than  the  P  velocity  in  the  mantle.  Other¬ 
wise  a  reversal  of  signs  takes  place  with  symme¬ 
try  about  normal  line. 


87 


3)  Larger  contrasts  of  velocities  between  the  cru3t 
and  the  mantle  result  in  larger  departures  of  the 
spectrum  of  the  apparent  angle  of  emergence  from 
normality.  Only  in  the  limiting  case  of  no  con¬ 
trast  of  velocities  between  crust  and  mantle  does 
the  apparent  angle  of  emergence  coincide  with  the 
normal  value  at  all  frequencies. 

4)  Rays  from  distant  teleseisms  incident  on  the  sur¬ 
face  at  less  than  25°  result  in  apparent  angles 
of  emergence  below  normal  when  very  large  period 
waves  are  considered.  Seismic  energy  radiated 
from  closer  sources  and  emerging  with  angles  of 
incidence  greater  than  3®°  has  an  apparent  angle 
of  emergency  which  slowly  increases  above  normal 
as  the  period  decreases. 

5 )  Maxima  of  the  apparent  angle  of  emergence  are 
controlled  by  the  minima  of  the  horizontal  com¬ 
ponent  transfer  function.  The  frequencies  of 
the  maxima  correspond  to  the  frequencies  at  which 
the  horizontal  component  of  motion  was  destroyed 
by  interference  patterns  of  the  primary  and  sec¬ 
ondary  waves.  The  energy  emerges  almost  vertically 
at  these  frequencies. 


6)  The  phase  of  the  apparent  angle  of  emergence  in¬ 
dicates  the  phase  shift  bctw*  \n  the  vertical  and 
horizontal  components  of  the  transfer  functions. 

At  zero  frequency  the  phase  is  180°.  In  a  par¬ 
ticle  motion  diagram  in  the  vertical  plane  this 
phase  relation  corresponds  to  the  linear  polariza¬ 
tion  proper  to  the  P  wave .  Following  Haskell 
(1962),  if  the  axes  be  chosen  as  indicated  in 
Figure  37 


Figure  37.  Components  of  ground  motion  for 

incident  P  wave. 

that  is,  horizontal  positive  in  the  direction  of 
advancing  wave  and  vertical  positive  downward,  the 
particle  will  vibrate  in  vertical  plane  and  with 
linear  motion.  As  the  frequency  of  the  wave  in¬ 
creases  the  phase  difference  between  the  two  com¬ 
ponent  i  changes  and  some  elliptical  polarization 


takes  place.  The  ellipticity  increases  with  the 

departure  of  the  phase  difference  from  the  normal 
o 

value  of  180  according  to  Figure  38 


Figure  38.  Particle  motion  diagram  of  P  wave. 

In  order  to  present  this  particle  motion  in  graphical 
form  a  FORTRAN  II  computer  program,  listed  in  Appendix  III, 
was  prepared.  The  output  is  presented  in  a  series  of 
graphs  (see  Figures  39  to  43).  The  particle  motion  is  in 
the  vertical  plane.  The  graphs  correspond  to  the  one-layer 
model  at  angles  of  incidence  of  35°>  and  45°  for  the 
contrasts  OU/  0£>  =  0.9  and  0.7. 

Not  only  does  the  ellipticity  of  the  particle  motion 
change  with  %  but  also  the  direction  of  the  major  axis  of 
the  ellipse.  The  ellipticity  of  the  ground  motion  is  great¬ 
er  for  large  contrasts  of  velocity  between  the  mantle  and 


Particle  motion  diagrams  of  ground  motion  at  different  values  or  the  parameter  f  for  Incident  P  wave  at  3**.  anp-le 
and  for  a  crustal  model  of  one  layer  where  j> *0.7 


Figure  40  :  Particle  motion  diorama  or  ground  motion  at  different  values  of  the  parameter  1  for  incident  P  wave  at  3^  an#rle  of  Incline 
and  for  a  crustal  model  of  one  layer  where  «t  •%/  ?■*  0.9 


article  motion  diagrams  of  ground  motion  at  different  values  of  the  parameter  ffor  Incident 
Lnd  for  crustal  model  of  one  layer  where  0.7  # 


the  ;rust.  Never the. less,  as  Indicated  by  Haskell  (1962), 
the  character  of  the  ground  motion  remains  generally  of 
the  P  wave  type. 

4.2  Two-layer  crustal  models 

The  introduction  of  an  intermediate  discontinuity  in 
the  crust  in  accordance  with  observations  in  many  differ¬ 
ent  parts  of  the  world  will  increase  the  number  of  reflec¬ 
tions  and  refractions  which  the  seismic  rays  undergo  in 
the  crust  of  the  earth.  Some  rays  will  be  refracted  at 
the  intermediate  discontinuity  and  so  their  amplitude  will 
be  affected.  The  time  lag  of  these  rays  will  also  be 
changed  and  the  periodicity  o.  the  transfer  functiona  will 
be  different.  Nevertheless,  observations  show  that  these 
discontinuities  have  relatively  small  velocity  contrasts; 
consequently  the  amplitudes  and  time  lags  Introduced  will 
be  small  and  the  deviation  from  the  one-layer  model  will 
not  be  as  great  as  might  first  be  expected. 

In  addition  to  the  above  rays  there  will  be  others 
due  to  reflections  at  the  intermediate  Interface.  In 
general  these  waves  will  have  reflection  coefficients  of 
small  amplitude.  Some  of  them  will  have  small  time  lags, 
as  for  example  the  reflections  of  P  at  the  interface  iden¬ 
tified  in  Figure  44  as  P2  Pi* 
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Figure  44.  Multiple  reflections  on  a  two-layer 

crustal  model. 

This  reflection  will  have  small  time  lags  as  compared  with 
the  direct  P,  and  consequently  will  introduce  a  long-period 
component  in  the  transfer  function  curve.  Other  reflections 
will  have  a  large  time  lag  and  their  contribution  will  in¬ 
troduce  short-period  sinusoidal  components  into  the  perio¬ 
dicity  cf  the  curves. 

For  the  two-layer  case,  many  combinations  of  the  layer 
thicknesses  and  velocities  ratios  are  possible.  Figures  45 
to  49  present  the  special  case  of  a  model  where  the  angle 
of  incidence  is  40°  and  c*-j/ =  0.9.  Each  family  of 
curves  corresponds  to  a  different  thickness  ratio  h-j/lv  . 

The  following  remarks  concerning  these  curves  may 
facilitate  their  application  to  observations: 


TRANSFER  FUNCTIONS  FOR  CRUSTAL  MODEL 


cruotaj  model* 


[•] 
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transfer 


TRANSFER  FUNCTIONS  FOR  CRUSTAL  MODEL 


In  the  crustal  model  h>  /  «-S ' 1  and 
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1)  In  general  the  two-layer  model  transfer  curves 

resemble  the  corresponding  one-layer  models.  This 
is  particularly  true  when  the  contrast  between 
the  P  velocities  in  the  crust  and  mantle  is  large, 
and  also  when  one  of  the  layers  is  very  thick 
with  respect  to  the  other.  For  example,  the 
curve  of  Figure  49  corresponding  to  ^  3/0^2  “  °*7 
and  h^/h,,  =  1/5  is  very  similar  to  the  one-layer 
case  =  0.7  presented  in  Figure  32. 

2)  The  transfer  function  of  models  with  layers  of 
the  same  thickness,  as  the  one  presented  in  Fig¬ 
ure  47  for  h^/h^  =  1/1 »  differ  considerably  from 
the  corresponding  one-layer  model. 

3)  The  large  amplitude  components  of  the  transfer 
•  functions  are  the  result  of  the  velocity  con¬ 
trast  betwc  1  the  mantle  and  the  crust. 

4)  The  influence  of  intermediate  layers  on  the  trans¬ 
fer  functions  increases  with  frequency.  This 
property  should  be  of  interest  vhen  the  spectrum 
of  short-period  observations  is  obtained. 

5)  The  periodicity  of  the  transfer  functions  within 

u  family  of  curves  is  similar  in  most  of  the  cases. 
Nevertheless  for  some  models  the  influence  of  the 
second  layer  is  so  small  due  to  low  contrast  of 
velocities  that  some  anomaly  in  the  periodicity 
of  the  curves  seems  to  be  present  (Figure  48). 
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6)  Changes  in  the  ratio  of  the  thicknesses  between 
the  two  layers  of  the  crust  Influences  the  trans¬ 
fer  functions.  This  influence  is  greater  if  the 
contrast  of  velocities  between  the  mantle  and  the 
crust  is  small.  Changes  of  thickness  ratios  af¬ 
fect  only  the  periodicity  of  the  transfer  func¬ 
tions.  This  is  true  because  these  changes  influ¬ 
ence  only  the  time  lags  of  the  different  arrivals 
at  the  surface  and  not  their  amplitudes. 
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5.  Spectral  Analysis  of  the  Seismic  Records 

As  previously  indicated,  the  tangent  of  the  apparent 
angle  of  emergency  can  be  obtained  from  the  seismograms  by 
dividing  the  spectrum  of  the  vertical  component  by  the 
spectrum  of  the  horizontal  component.  This  operation  will 
eliminate  from  the  spectra  the  influence  of  the  frequency 
content  of  the  source. 

The  theoretical  justification  for  this  property  of 
the  apparent  angle  of  emergence  is  given  by  Hannon  (1964) 
and  by  Phinney  (1965). 

To  be  completely  valid  the  argument  assumes  that  for 
a  seismic  ray  the  frequency  content  at  the  base  of  the 
crust  is  exactly  the  same  for  both  horizontal  and  vertical 
components.  In  other  words,  before  the  mantle-crust  dis¬ 
continuity,  no  other  significant  discontinuity  was  en¬ 
countered  by  the  ray.  This  is  only  partially  true  since 
the  upper  mantle  and  mantle  are  not  exactly  homogeneous. 
Even  so,  it  may  reasonably  be  assumed  that  deep  in  the 
mantle  velocity  gradients  exist  and  that  changes  in  the 
velocity  of  P  waves  are  gradual;  no  reflections  from  dis¬ 
continuities  can  be  expected  for  this  part  of  the  ray-path 
and  the  transfer  function  for  this  portion  of  the  medium 
will  be  unity,  according  to  the  theoretical  considerations 
of  section  3.  Should  discontinuities  exist  in  the  upper 


mantle  the  velocity  contrasts  across  these  discontinuities 
may  reasonably  be  presumed  to  be  relatively  small;  quali¬ 
tatively,  in  view  of  the  previous  discussion,  the  effect  of 
these  discontinuities  would  be  to  introduce  a  small  oscil¬ 
latory  component  into  the  form  of  the  transfer  function. 
Because  the  velocity  contrasts  are  small  the  amplitude  of 
this  component  would  be  small.  Because  the  time  lags  so 
introduced  would  be  large,  the  periodicity  of  this  compo¬ 
nent  would  be  small.  The  net  result,  in  general,  would  be 
a  small  rapidly  oscillating  character  introduced  into  the 
form  of  the  transfer  functions  riding  on  top  of  the  main 
features  of  these  functions  due  to  discontinuities  of  the 
crust. 

An  investigation  of  the  upper  mantle  structure  could 
be  attempted  based  on  the  observations  we  have  just  made 
in  the  present  study;  however,  we  do  not  intend  to  follow 
this  procedure.  In  fact,  we  shall  smooth  the  spectra  to 
be  obtained  from  seismograms,  and  one  reason  for  so  doing 
will  be  to  exclude  any  oscillatory  character  of  the  trans¬ 
fer  function  of  short  periodicity  due,  possibly  to  the 
influence  of  the  upper  mantle . 

5.1  Fourier  analysis  of  longitudinal  seismic  waves. 

In  the  spectral  analysis  of  the  seismograms  several 
restrictions  are  imposed  by  the  very  nature  of  spectral 
analysis  and  by  the  transient  aspect  of  the  seismic  body 
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The  finite  length  of  the  seismic  record  limits  the 
frequency  resolution  of  the  spectrum.  This  resolution  is 
a  function  of  the  total  length  of  the  seismogram  used  to 
estimate  the  spectrum.  The  resolution  is  given  by  Blackman 
and  Tukey  (1958)  by  the  relation 

(resolution  in  cps)  =  - i - 

(total  length  in  seconds) 

This  notion  of  resolution  is  intrinsically  connected 
with  the  concept  of  a  data  window  and  with  smoothing  of  the 
spectrum.  In  practice,  when  a  finite  length  of  record  is 
taken,  if  no  shaping  or  tapering  of  the  record  is  intro¬ 
duced  the  record  is  actually  passed  through  a  rectangular 
window  of  the  form 

Do  (t)  =  1,  for  t  C  Tm 
=  0,  for  t  ^  Tm 

where  Tm  is  the  total  length  of  the  record  available. 


h 

- 1 - 

i 

i 

- 1 - — p. 

-T*  0  Tr»  t. 

Figure  50.  Rectangular  time  window. 

When  this  record  is  Fourier  analyzed  the  result  is  the  con 
volution  of  the  true  spectrum  P  (^  )  of  the  signal  and  the 
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Fourier  transform  of  the  rectangular  window  Qq  (  |  ) .  This 
spectrum  is  only  an  estimation  of  the  true  spectrum  and 
may  be  called  the  "average  spectrum" 


Consequently,  the  value  estimated  for  the  spectrum  at  the 
frequency,  f^,  is  the  average  value  at  this  frequency  and 
the  neighboring  frequencies  with  a  weight  proportional  to 
Qq  (f^-f).  The  Fourier  transform  of  the  rectangular  func 
tion  is 


2  Tm 


sm  (2 I**) 

zn  J  T* 


(5-3) 


and  its  graph  is  given  in  Figure  51. 


Figure  51.  Fourier  transform  of  the  rectangular 

window 
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This  figure  illustrates  the  fact  that  if  no  special  data 
window  is  used  in  the  analysis  of  the  records,  in  reality 
a  rectangular  window  is  used  and  the  spectrum  will  be  a 
poor  estimation  of  the  real  spectrum.  This  follows  since 
the  spectral  window  does  not  behave  in  the  manner  we  might 
desire;  the  estimation  at  each  frequency  is  affected  by 
the  spectral  amplitude  values  of  distant  frequencies. 

Several  data  windows  have  been  devised  to  remedy,  at 
least  in  part,  this  difficulty.  For  example,  Kirasawa  and 
Stauder  (1964)  use  a  data  window  of  exponential  type  pre¬ 
sented  by  Kasahara  (1957).  Phinney  (1964)  uses  for  his  lag 
window  a  Gaussian  function.  For  the  present  analysis  we 
have  selected  the  use  of  the  taper-pair  window  given  by 
Blackman  and  Tukey  (1956)  and  designated  by  them  as  a 
"hamming"  window.  Its  form  is 


lfc|  <  X 

M  >  T*  (5'4) 

°>(j)  *  t5  *'4 (£) «-  «■  «  (J+  If-  x 

where  Qq  iS  given  by  (5-3). 

The  graph  of  the  pair  is  given  in  Figure  52.  This  figure 
shows  that  for  a  sufficiently  large  value  of  Tm  the  estima¬ 
tion  of  the  frequency  is  restricted  to  a  narrow  frequency 


0  ft)  -  O.Sr^t  +  0.  4  6  cos  JL 

~T 

-  (9 

and  its  Fourier  Transform  is: 


band  without  influence  of  the  side  lobes . 


Figure  52.  Hamming  data  window  and  its  Fourier 

Transform 

In  addition  to  the  smoothing  aspect  of  the  data  window 
its  use  serves  another  purpose.  Tine  data  window  can  be  con¬ 
sidered  as  a  weighting  function  which  emphasizes  the  first 
part  of  the  P  wave  and  minimizes  the  influence  that  later 
arrivals  could  have  on  the  spectrum.  This  is  very  conven¬ 
ient  in  the  analysis  of  the  seismic  body  waves  since  the 
first  arrival  of  the  P  wave  is  followed  at  short  time  Inter¬ 
vals  by  other  phases  such  as  pP  or  PcP  not  related  to  the 
direct  P  wave.  In  practice  the  selection  of  the  total  length 
of  the  record  to  be  used  and  the  value  of  the  window  for 
which  the  weight  should  be  zero  is  very  critical.  If  too 
short  a  time  is  selected  some  important  reflections  will 
not  be  included  and  the  transfer  function  will  be  incomplete. 
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If  too  long  an  interval  is  selected  other  phases  as  pP  and 
PcP  could  be  included  and  the  transfer  function  obtained 
will  no  longer  be  representative  of  the  crustal  structure 
of  the  station  under  study. 

Proper  selection  of  the  time  interval  of  record  to  be 
used  is  a  function  of  the  magnitude,  depth,  and  epicentral 
distance  of  the  earthquake.  Some  estimation  of  the  mini¬ 
mum  time  required  may  be  obtained  from  the  calculations  of 
Ben-Menahem,  et  al.  (1964),  who  computed  the  theoretical 
record  corresponding  to  different  structures  and  a  conven¬ 
tional  source  function.  The  computed  seismograms  include 
the  reverberations  from  the  layers  of  the  cri'st  and  show 
that  at  times  greater  than  40  seconds  after  the  first 
arrival  there  are  no  longer  substantial  contributions  to 
the  record  from  the  direct  P. 

Another  limitation  of  spectral  analysis  is  introduced 
by  the  folding  frequency .  To  avoid  this  folding  of  the 
spectrum  it  is  convenient  to  use  linear  digital  filters  with 
a  low  pass-band.  This  may  be  done  using  one  of  any  number 
of  filters.  We  shall  here  use  a  digital  linear  smoothing 
filter  with  no  phase  shift  presented  by  Hamming  (1962). 

This  filter  replaces  the  original  value  of  the  record  yk(t) 
by  a  new  one  given  by 

^  *  JM  +  2  M  f  Mu 

_  k _ j  n*  1  J  _ J*-4 1  Jn4<, 

l  2. 


(5-5) 


1.1,1 


The  ability  of  thio  filter  to  reject  the  upper  half 
of  the  spectrum  is  very  good  and  its  transfer  function  is 


where  T  is  the  period  of  the  wave.  This  filter  puts  zeros 
in  the  spectrum  at  1/4,  1/3  and  1/2  of  the  folding  frequency 
as  is  shown  in  Figure  53. 


Figure  53.  Frequency  response  or  digital  filter. 

The  transient  character  of  the  P  wave  does  not  allow 
us  to  treat  it  as  a  stationary  signal.  Consequently,  in- 
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stead  of  calculating  the  autocorrelation  of  the  signal  and 
the.'  taking  the  Fourier  transform  to  find  the  power  spec¬ 
trum,  the  records  are  directly  Fourier  integrated*  The 
numerical  integration  is  performed  using  Filon's  method  as 
described  by  Hamming  {1962). 

In  order  to  combine  all  these  operations  In  a  flexible 
FORTRAN  program  a  main  program  was  written  which  includes 
e.  series  of  subroutines  with  the  different  operations. 

This  program,  called  Seismic  Analyzer,  is  listed  in 
Appendix  IV.  The  output  of  this  program  is  the  series  of 
graphs  used  in  the  next  section  for  crustal  determinations. 
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6.  Crustal  Parameters  of  the  Central  United  States 


In  this  section,  as  an  application  of  the  theory  and 
methods  of  analysis  developed  previously,  we  shall  eva3uate 
the  thickness  of  the  c"ust  and  the  average  crustal  P  wave 
velocity  for  the  region  of  the  central  United  States.  The 
data  for  this  evaluation  will  consist  of  the  spectra  of  tne 
apparent  angle  of  emergence  obtained  from  seismograms  re¬ 
corded  at  stations  of  the  Saint  Louis  University  network. 

The  interpretation  itself  will  progress  in  two  steps. 
First,  the  observational  spec era  will  be  compared  to  theo¬ 
retical  spectra  for  single-layered  models  of  the  crust.  In 
this  way  a  first  approximation  to  the  crustal  parameters 
may  be  obtained.  Second,  two- layered  models  snould  be  used 
in  an  attempt  of  a  more  precise  interpretation. 


6 . 1  Application  of  the  one-layer  model  master  curves. 

For  the  one-layer  models  the  dimensionless  parameter  Y 
becomes 


Since  the  families  of  master  curves  were  calculated  for 
constant  values  of  f^  and  e1  it  is  possible  to  divide  Y by 
tne  expression  in  the  parenthesis.  Then 


(6-1) 
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This  simplification  of  represents  a  translation  of  the 
curves  along  the  ^  axis  when  they  are  plotted  on  loga¬ 
rithmic  paper.  A  sot  master  curves  with  this  transla¬ 
tion  is  presented  in  Figures  54  to  62*.  In  these  figures 
only  the  tangent  of  the  apparent  angle  of  emergence  is 
given.  Each  family  of  curves  corresponds  to  an  angle  of 
incidence  of  the  ray  in  the  crustal  layer.  Only  three  mem¬ 
bers  of  the  family  are  presented  for  the  sake  of  simplicity 
and  clarity.  Intermediate  values  of  the  ratio  may  be 
easily  interpolated. 

The  use  of  these  master  curves  may  be  summarized  by 
the  following  procedural  steps: 

l)  Evaluate  the  tangent  of  ‘one  apparent  angle  of  emer¬ 
gence  versus  frequency  at  a  given  station  using  the  method 
and  program  presented  in  the  preceding  chapter. 

The  selection  of  the  earthquakes  for  this  type  of 
analysis  is  very  critical.  The  main  criterion  is  that  at 
least  40  seconds  of  the  seismogram  correspond  to  the  direct 
P  recording  free  of  interference  of  other  phases  which  have 
their  origin  outside  cf  the  crust  below  the  station.  In 
practice  the  phases  to  be  avoided  are  pP  ana  PcP,  since 
these  are  the  waves  which  arrive  at  a  short  time  interval 
after  the  direct  P  phase.  Records  which  satisfy  this  con¬ 
dition  can  be  obtained  at  epicentral  distances  smaller 


*To  facilitate  the  use  of  these  curves,  they  are  included 
in  a  pocket  at  the  back  cover. 
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than  55°  and  for  earthquakes  with  a  depth  of  focus  greater 
than  150  km.  The  magnitude  of  the  shocks  should  be  such 
that  the  signal-to-noise  ratio  of  the  long-period  records 
is  sufficiently  high. 

2)  Plot  the  curve  of  the  tangent  of  the  apparent  angle 
of  emergence  versus  frequency  on  logarithmic  tracing  paper 
of  the  same  scale  as  that  of  the  master  curves. 

3)  Rrom  the  available  information  of  epicentral  dis¬ 
tance  and  depth  of  the  focus  estimate  roughly  the  angle 

of  incidence  of  the  ray  at  the  surface  of  the  earth.  This 
evaluation  is  only  a  guiding  criterion  in  the  selection  of 
the  correct  family  of  maste  '  curves  to  be  used. 

4)  Match  the  observed  and  theoretical  curves.  To  ob¬ 
tain  this  matching  the  observed  curves  are  laid  over  the 
theoretical  curves  and  then  moved  to  the  right  and  to  the 
left,  up  and  down,  but  without  rotation  of  the  axis  until 
the  best  match  is  achieved.  In  the  matching  special  atten¬ 
tion  should  be  given  to  the  central  portions  of  the  ob¬ 
served  curves  since  these  correspond  to  the  central  values 
of  the  response  curve  of  the  seismometer  lor  which  the  ob¬ 
servations  are  more  accurate. 

5)  The  value  of  the  velocity  contrast  may  be  obtained 
by  the  interpolation  between  the  members  of  the  families 
of  curves.  If  the  velocity  of  the  mantle  in  the  region 

is  known  from  other  sources  or  assumed  in  some  way,  the 
average  velocity  of  the  P  waves  in  the  crust  0(1  rnay  be 


obtained  from  the  ratio  ol^/ol  2» 

6)  Read  the  value  of  $  corresponding  to  the  0.1  cps 
frequency  value  of  the  observations.  If  we  call  this 
value  ,  then: 


(6-2) 


and  the  thickness  of  the  crust  will  be: 

4  =  10.  1 Tv 

A  graphic  example  of  the  method  is  presented  in  Figure  63. 


0.1  0.2  0.3  0.4  0.6  1.0  2.0 

Figure  63.  Matching  of  theoretical  and  observed 

curves . 
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In  Figure  63 

=  0.82  oLx 

1 

X,  =  0.61 

=  0.6l  x  10  x  0.82  =  41.5 

Go  assuming  a  value  of  the  P  velocity  In  the  mantle  of  8.2 
km/sec.,  the  average  P  velocity,  =  6.806  km/sec.  and 
the  total  thickness  of  the  crust  h  =  41.5  km. 

6.2  Multi-layered  cas-s. 

Application  of  the  master  curves  for  the  case  of  two 
and  more  layer  models  may  be  done  in  a  similar  way. 

In  the  case  of  two  layers  it  is  possible  to  transform 
the  dimensionless  parameter  as  follows: 


e.-f  >/j  jiUJ,) 


* - — (jl«  ) 

(6-31 


otA 


Since  for  each  family  of  master  curves  h^/h^  and  ot0/o( ^ 

remain  constant,  it  Is  possible  t^  divide  by  the  quantify 

! 

in  brackets  and  present  the  curves  in  terms  of  X  =  h 
In  this  way  the  thickness  h-,  and  the  velocity  of  the 
first  layer  may  be  obtained  directly  from  the  matching  of 
the  curves.  We  then  calculate  h^  and  for  the  second 
layer  from  the  ratios  h^/hg  and 

A  set  of  90  master  curves  for  two-layer  models  have 


1* 


been  calculated.  This  set  Includes  three  velocity  con¬ 
trasts  in  the  crust  with  values  for  01  0(2-0.90,  O.95 


and  1.0b.  The  third  case  corresponds  to  a  low  velocity 
layer  't  the  bottom  or  the  crust.  The  ratio  of  tne  layers 
is  h1/h2  =  1/5,  1/2,  1/1,  2/1,  5/1. 

6 . 3  Application  to  the  crustal  structure  in  Central 
United  States. 

The  master  curve  matching  technique  just  described 
was  applied  to  long-period  seismograms  obtained  at  the 
seismic  network  of  stations  of  Saint  Louis  University  in 
the  central  United  States.  The  list  of  stations  and  their 
coordinates  is  given  in  Table  1 . 

Table  1 .  STATIONS  OF  SAINT  LOUIS  UNIVERSITY  NETWORK 
Geographic  Latitude  Longitude 

Bloomington  39°  11 «  20"N  86°  30'  1S"W 

Dubuque  42°  30'  24"N  90°  4l«  00,!W 

Florissant  38°  48'  06"N  90°  22'  12" W 

Manhattan  39°  11'  59"N  96°  34'  50"W 

Rolla  37°  55 !  04"N  91°  52'  08"W 

Saint  Louie  38°  38'  10"N  90°  14'  10"W 

The  response  characteristics  of  the  long-period  seismograph 
used  in  this  study  were  calculated  by  Nuttli  and  McEvilly 
(1961)  using  the  impedance  bridge  method  described  by  Will- 
more  . 

According  to  the  criteria  previously  indicated,  three 
earthquakes  of  intermediate  and  deep  focus  were  selected. 
The  three  components  of  the  long-period  instruments  were 
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used  to  evaluate  the  spectrum  of  the  angle  of  emergence 
according  to  the  method  indicated  in  section  9  of  this 
study.  The  parameter  of  the  three  selected  earthquakes 
according  to  the  U.S.  Coast  and  Geodetic  Survey  Bulletins 
are: 


Date  Origin  time  Epicenter 
Year  Mo.  Day  Hr  .min. sec .  Lat.  Long.’ 


Depth  Magn.  Site 

PT 


1961  Dec.  20  13-25-34.4  4.6n 
1963  Aug.  15  17-25-05.9  13.8S 
1963  Nov.  9  21-15-30.4  9.0S 


75. 6w 

176 

6  3/4  Colombia 

69. 3W 

543 

7  3/4  Peru 

7.1W 

600 

6  3/4-7  Brazil 

The  three  long-period  components  were  digitized  at  a 
rate  of  one  sample  every  O.566  seconds,  which  corresponds 
to  a  folding  frequency  of  1  cps.  The  length  of  the  time 
window  in  each  case  was  determined  from  the  epicentral  dis¬ 
tance  and  depth  of  the  focus . 

The  signal-to-noise  ratio  of  some  seismograms  was  such 
that  the  analysis  could  be  performed  without  the  filtering 
and  detrending  steps . 

The  result  of  the  final  analysis  is  presented  in 
table  2  and  in  Figures  64  to  75.  In  these  figures  the  orig¬ 
inal  records  are  presented  together  with  the  appearance  of  the 
phase  after  the  filtering  and  weighting  operations  have 
been  performed.  Pt  the  left  side  of  the  graph  the  tangent 
of  the  apparent  angle  of  emergence  versus  frequency  is  given 
together  with  the  ohase  difference  between  vertical  and  hori¬ 
zontal  components  of  the  Fourier  analysis. 


APPARENT  ANGLE  OF  EMERGENCE  VS.  FREQUENCY 
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FREQUENCY  (CPS)  Figure  64  :  Crustal  solution  at  Blootaliigton  for  earthquake  of  August  IS.  19t 


APPARENT  ANGLE  OF  EMERGENCE  VS.  FREQUENCY 


FREQUENCY  ICPW)  Pl«ure$5  :  Cruatal  ablution  at  Bloowtng.on  for  the  earthquake  of  November  9,  19f 3  • 


APPARENT  ANGLE  OF  EMERGENCE  VS.  FREQUENCY 


FPEOUENC Y  ICPS)  Plgure  Vr  :  Crustal  solution  at  PlorlBsant  for  earthquake 


APPARENT  ANGLE  OF  EMERGENCE  VS.  FREQUENCY 


Crustal  solution  at  Manhattan  Tor  earthquake 


APPARENT  ANGLE  OF  EMERGENCE  VS.  FREQUENCY 


Figure 71  :  Crustal  solution  at  3alnt  Louis,  for  earthquake  of  Novemter  ?,  1  . 


APPARENT  ANGLE  OF  EMERGENCE  VS.  FREQUENCY 
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FREOUENCYfCPT)  Pl&ura  ?3  :  Cruatal  solution  at  Rolla  for  earthquake  of  Auauat  15.19t} 
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1 


Flgupa  lit:  Cruatal  solution  at  Holla  for  earthquake  of  November  9,  19t3. 
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The  interpretation  of  the  observed  curves  given 
in  Table  2  gives  very  consistent  values  for  the  one-layer 
models.  Considering  the  area  covered  by  the  stations  as 
a  unit,  an  average  P  velocity  in  the  crust  of  6.6  Km/sec 
is  obtained.  The  average  thickness  of  the  crust  is  42  Kjn., 
The  above  results  are  in  good  agreement  with  values 
of  the  crust  obtained  by  previous  investigators  using  dif¬ 
ferent  techniques.  A  summary  of  these  studies  is  presented 
by  McEvilly  (1964)  and  is  reproduced  here: 


SUMMMARY  OF  CRUSTAL  STUDIES  IN  THE  CENTRAL  U.S. 

Crustal 


Reference 

Observations 

Pg 

Pn 

Thickness 

Walter  (1940) 

Earthquakes  near 
St.  Louis 

6.03 

7.73 

39 

Stechschulte  & 
Birkenhauer 

Ohio  earth¬ 
quakes 

8.25 

(1948) 

Lehman  (1955) 

N.E.U.S.  &  Okla¬ 
homa  earthquakes 

8.16 

Nutt 11  (1956) 

Tennessee  earth¬ 
quake 

6.1 

8.19 

Steinhart  & 
Meyer  (1961) 

Ark. -Mo.  Refrac¬ 
tion  profile 

5.18 

8.16 

4l 

Stauder  & 
Bollinger 

(19*3) 

McCamy  &  Meyer 
(1963) 

S.E.  Missouri 
Earthquake 

6.41 

8.20 

— 

Mo. -Ark.  Refrac¬ 
tion  profile 

6.2 

8.1 

45-50 

McEvilly  & 
Stauder 
(1963) 

Body  and  Ray¬ 
leigh  waves  from 
strip-mine  ex¬ 
plosions,  Ill. 

6.1 

Stewart  & 
McEvilly 

USGS  Mo.  Refrae- 
tion  data — pre¬ 
liminary  inter¬ 
pretation 

6.1 

8.15 

36-40 

3  3  4 


McEvilly,  using  phase  and  group  velocities  of  surface 
waves,  obtains  a  model  of  38  Km  thickness  and  an  average  P 
velocity  of  6.5  cm/sec. 

The  interpretation  of  the  observed  curves  in  terms 
of  two  layer  crustal  models  is  subject  to  more  ambigu¬ 
ities  and  different  models  may  be  selected  to  match  the 
observations.  This,  of  course,  indicates  that  the  inter¬ 
pretation  given  in  Table  2  for  the  two-layer  case  is  not 
unique.  Nevertheless  all  the  different  interpretations 
give  a  very  constant  value  for  the  total  thickness  of  the 
crust.  The  average  value  of  this  thickness  for  the  Central 
United  States  using  the  two-layer  models  is  42  Km. 

It  is  expected  that  the  use  of  short-period  records 
and  master  curves  corresponding  to  the  0.1  and  1.0  cps 
frequency  band  will  provide  more  detailed  information  on 
the  layering  of  the  crust  and  a  better  interpretation  of 
observational  data  in  terms  of  two  or  more  layer  models. 


7.  Summary  and  Conclusions 
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1)  The  use  of  the  aimens ionless  parameter  allows  the 
systematization  of  the  effects  of  layered  systems  on 
longitudinal  seismic  waves  into  families  of  transfer 
function  curves. 

2)  The  periodicity  of  the  transfer  function  curves  is  de¬ 
termined  by  the  relative  time  lags  of  the  arrivals  at 
the  surface  of  the  reflections  and  refractions  within 
the  layered  system. 

3)  The  amplitude  of  the  transfer  function  curves  is  de¬ 
termined  by  the  contrasts  of  the  P  wave  velocities  at 
the  interfaces  of  the  layered  system. 

4)  A  set  of  normalized  master  curves  of  the  transfer  func 
tions  for  one- and  two-iayer  models  and  for  several 
angles  of  incidence  has  been  prepared. 

5)  The  use  of  deep  and  intermediate  depth  earthquakes  at 
convenient  epic.entral  distance  allows  the  calculation 
of  the  spectrum  of  the  apparent  angle  of  emergence 
independently  of  the  conditions  at  the  scurce.  Tni0 
spectrum  depends  only  on  the  structure  Immediately 
below  the  receiving  station. 

6)  In  most  of  the  cases  the  matching  of  master  curves 
and  observations  allows  the  calculation  of  the  thick¬ 
ness  and  average  P  velocity  of  the  crust  with  less 
than  +  5$  incertitude. 


7)  Estimation  of  different  layers  within  the  crust  from 
the  spectrum  of  long-period  records  is  neither  precise 
nor  unique.  It  is  suspected  that  spectra  of  short- 
period  records  will  give  more  detailed  information  of 
the  structure  within  the  crust. 

8)  Application  of  the  method  to  the  region  of  central 
United  States  occupied  by  the  Saint  Louis  University 
network  of  seismic  stations  gave  an  average  crustal 
thickness  of  42  Km  find  average  crustal  velocity  of  F 
waves  of  6.6  Km/sec. 

9)  The  U3e  of  the  dimensionless  parameter  $  should  facil¬ 
itate  the  calculation  of  the  first  partial  derivatives 
of  the  transfer  functions  wich  respect  to  the  parameters 
of  the  layers.  Depending  on  the  behavior  of  these  par¬ 
tial  derivatives  the  inversion  problem  for  the  spectrum 
of  F  body  waves  may  be  attempted.  This  inversion  prob¬ 
lem  r hould  find  the  best  fit,  in  a  least  squares  sense, 
between  observations  and  theoretical  curves. 
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**  FERNANDEZ  6  PHI  8  TRANSFE 

C  THIS  PROGRAM  IS  DESIGNED  TO  CALCULATE  THE  TRANSFER 

C  FUNCTIONS  OF  A  LAYERED  SYSTEM  FOR  THE  VERTICAL  AND 

C  HORIZONTAL  COMPONENTS  OF  dilatational  WAVES  AND 

C  THEIR  RATIO,  THIS  RATIO  IS  THE  TANGENT  OF  THE 
C  APPARENT  ANGLE  OF  EMERGENCE. 

C  THE  TRANSFER  FUNCTIONS  MAY  BE  CALCULATED  IN  TERMS 

C  Or  FREQUENCY  OR  IN  TERMS  OF  A  DIMENSIONLESS 

C  parameter  gamma. 

C  THE  PROGRAM  WILL  PLOT  ALSO  THE  MODULUS  AND  THE  PHASE 
C  OF  THE  TRANSFER  FUNCTIONS  IN  ORDINARY  OR  LOGARITHMIC 
C  SCALE.  THE  VALUES  OF  THE  TRANSFER  FUNCTIONS  OF  THE 
C  PLOT  ARE  NORMALIZED  TAKING  AS  UNITY  THE  CORRESPONDING 

C  COMPONENT  OF  MOTION  AT  THE  BOTTOM  OF  THE  LAYER  SYSTEM. 

C  THE  APPARENT  SURFACE  VELOCITY  MUST  BE  GREATER  THAN 
C  THE  P  VELOCITY  IN  THE  TOP  LAYER  AND  THE  S  VELOCITY 
C  IN  ANY  LAYER. 

C  THE  ANGLE  OF  INCIDENCE  OF  THE  INPUT  CORRESPONDS  TO 

C  THE  ANGLE  OF  THE  RAY  AT  THE  SURFACE. 

DIMENSION  D(50),A(50)#B(50),RH0(50),DGAM(50) 

C  READ  THE  SCALE  FACTORS  AND  CONSTANTS  FOR  THE  PLOT  OF 
C  THE  FUNCTIONS. 

READ  1 92#XS,YS,YPS,X02 
READ  1 92,Y01 ,Y02,Y03,Y04 
192  F0RMAT(4F10.0) 

1  030  PUNCH  1 

1  F0RMAT(1X,!8HTRANSFER  FUNCTIONS) 

C  UNDER  STUDY 

C  READ  THE  IDENTIFICATION  OF  THE  CRUSTAL  MODEL 

READ  1101 
PUNCH  1101 

110!  F0RMAT(1X,20H  ) 

C  NOL-NUMBER  OF  LAYER5+1 ,  AT-1  IF  GAMMA  DES IRED,PL0TR-1 
C  IF  PLOT  OF  RATIO  ONLY  DESIRED,  PLOTO-1  IF  PLOT  NOT 
C  IN  LOGARITHMIC  SCALE  DESIRED 

READ  2.N0L,AT,PL0TR,PL0T0 

2  FORMAT* I 3.3F3.0) 

PUNCH  2.N0L,AT,PL0TR,PL0T0 
PUNCH  707 

707  F0RMAT(10HTHICK.(KM),1X,14HP  VEL. (KM/SEC), 1X.14HS 
1  VEL, (KM/SEC),  IX. 1 2HDEN.(GM/CMC),1  X,5HLANDA,/) 

C  READ  THE  LAYER  PARAMETERS.  D-THICKNESS,  A»P  VELOCITY, 

C  B-S  V£LOCITY,RHO-  DENSITY 

DO  3  1*1  ,NOL 

READ  500,D(l),A(!)fB(l),RH0(l) 
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500 

3 

501 


C 

C 

C 

7001 

4 

C 

C 

c 

c 

51  04 
831 

5004 


742 

741 


1346 


237 


286 

287 

285 


F0RMAT(1X,F9.4.3F10.4) 

FLAND( I ) *RH0( I )*A( I )*A( I ) /3. 

PUNCH  501  ,D(  I  ),A( !  ),B(  l),RH0(  I  ),FlAND(  I ) 

F0RMAT(5F1  0.4) 

u*N0L 


N1  =N0L“1 
N2*N0L-2 

READ  ANGlIelNITlAL  ANGLE  OF  INCIDENCE  AT  THE  SURFACE 
DANG ! =  INCREMENT  OF  THE  ANGLE  OF  INCIDENCE  , 
NOANG=NUMBER  OF  ANGLES  OF  INCIDENCE  TO  BE  CALCULATED 
READ  4.ANGI1 ,DANGI,NOANG 
FQRMAT(2F10.5,I3) 

GMIN*MINIMUM  VALUE  OF  GAMMA  DESIRED,  GINC-  INCREMENT 
OF  GAMMA, N OF*  NUMBER  OF  VALUES  OF  GAMMA  DESIRED. 

IF  THE  CALCULATIONS  ARE  IN  TERMS  OF  FREQUENCY  „ 

LET  GMIN*  MINIMUM  FREQUENCY  DESIRED,  ETC. 

READ  4,GM1N,GINC,N0F 

ANG I P*ANG 1 1  -DANG  I 

DO  311  IA=1,N0ANG 

GAF*£MIN-GINC 

ANG  I  P*ANG  I P+DANG  I 

S IN  II *SINF(ANG IP/57.29578) 

SININ=A(N0L)*SINI1/A(1  ) 

C=A(NOL)/S ININ 

S!NE*(1  .**2.*(S  ININ**2)/3.) 

CCSF*SQRTF(1 .-S IN£**2) 

TANE=*SINE/COSE 
£M*ATANF ( TAN  E)*57. 29578 
AIM*90.-EM 
PUNCH  742 

FORMAT ( 9HAP . VEL OC . , 2X ,  1  OHTAN . EMERG . ,  1  4HANG .  INC. 

1  30T. ,1 5HANG.  INC.  SURF.) 

PUNCH  741  ,C,TANE,A IMjAHG IP 

FORMAT (4Fi 2.3) 

DO  1346  M=»1,N0L 
C0VA*C/A(M) 

COVB*C/B(M) 

DGAM(MH2./(C0VB**2  ) 

DGAM1 (M)«0GAM(M)-1 . 

DIFLA(M)-  SQRTF(ABSF(C0VA**2-  1.)) 

DRB(M)  «SQRTF(A6SF(C0VB**2-1 .)) 

DH(M)  «RHG(M)*C*C 
TH»0.  • 

APTH*C. 

PUNCH  TAPE  237 


F0RMAK1H*) 

IF  (AT 3285, 286, 286 
DO  28?  L«1  ,N1 
TH-TH+D(L)/A(L) 

APTH-APTH4D(L)*(SQRTF(1.-(A(L)*S!NI1/a(1  ))**2)  + 
1  1.732005*SQRTF(1 .-(B(L)*SINI1/a(1  ))**2))/a(L) 
CONTINUE 
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PUNCH  499 

lm  FORMAT (1  OX, 28HC0NSTANTS  OF  PROPORTIONALITY,/) 

PUNCH  708 

708  F0RMAT(5HLAYER,1X,11HTHICK.RATI0,1X,7HP  RAT  10,1  X, 
WHS  RATIO, 2X, 1  OHDEN.  RATI0,2X,7HL  RATIO,/) 

DO  502  I  =1 , NOL 
PK=D{ i)/TH 
QK«A( I ) /A (NOL) 

VK=B{ : )/B(NOL) 

RK=RHG  (  i )  /  RHO  ( NOL ) 

SK»FLAND( !)/FLAND(NOL) 

502  PUNCH  503, I ,PK,QK,VK,RK,SK 

503  FORMAT (1  X, 1 2,5X,F9.4,4F1 0.4,/) 

PUNCH  1620 

1620  FORMATS X,6hPERIOD,3X,4HFR£Q,3X,4HHF/A,3X,2HTW,5X, 
l3HPHV,',6x,2HTU,5X,3HPHU,3X,5HTW/TU,lXt4HNORM,2X, 
15HPHASE) 

DO  310  IFR-1.N0F 
GAF=*GAF+G  INC 
IF(AT)  290,291,292 

290  ■  FR£Q=GAF 

GO  TO  293 

291  FREQ=GAF*A(NGL)/TH 
GO  TO  293 

292  FR£Q=GAF/APTH 

293  CONTINUE 
PER-1 ./FREQ 

WVN 0=6. 2831 853*FREQ/C 

Al  1  =1 

Al  2=0 

A21  =0 

A  22-1 

A31"*0 

A32=0 

A4l  =0 

A42=0 

C  COMPUTE  ELEMENTS  OF  A  MATRIX  FOR  REMAINING  LAYERS 
DO  1345  M«1,N1 
GAM=UGAM(M) 

GAMM1  =DGAM1  (M) 

ra=dra(m) 

rb=drb(m) 

H=DH{M) 

P=WVNO*D(M)*RA 
Q=WVNO*D ( M) *RB 
SINP=SINF(P) 

W=S INP/RA 
X=RA*S  I NP 
C0SP»C0SF(P) 

S INQ=S INF(Q) 

Y«S!NQ/RB 

Z«RB*SINQ 

COSQ»COSF(Q) 

RHOH=RHO(M) 

BM“B(M) 
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DM=D(K) 

Bl 1 «GAM*C0SP-GAMM1  *COSQ 
31  2«GAMM1 *W+GAM*Z 
B1 3““(C0SP-C0SQ)/H 
614=  (W+Z)/H 
621  =  GAM*X+GAMM1  *Y 
B22*=~GAMM1  *COS  P+GAM*COS  Q 
B23  =“(X+Y)/H 
324  =B1  3 

83 i  bh*gam*gammi*(cosp-cosq) 

B32  =H*(GAMM1  *GAMM1  *W+GAM*GAM*Z) 

333  =B22 
B31'  =B1  2 

3 4 .  «-*H*  (  G AM*G  AH*X  +G  AMM1  *GAMM1  *Y ) 

R42  =B31 
B43  =B21 
B44  =B11 

EAl  1  =B1 1  *A1 1  +  B1  2*A21  +  B1  3*A31  +  91  4*a41 
EA1  2  «  Bl  1  WA1  2  +  Bl  2*A22  +  Bl  3*A'*>2  +  B1  4*a42 

EA21  «  B21*A11  +  322*A21  +  B23*A31  +  B24*A4l 

EA22  «=  B21*A12  +  B22*A>2  +  B23*A32  +  B24*a42 

EA31  -  B31 *A1 1  +  L32*A21  +  B33*A31  +  B34*a41 

EA32  «  B31*A12  +  B32*a22  i-  B33*A32  +  B34*a42 

EA41  «  B41 *Al 1  +  B42*a21  +  B43*A31  +  B44*a41 

EA42  »  B41 *A1 2  +  B42*A22  +  B43*A32  +  B44*a42 

Al  1  «  EAl 1 
Al  2  »  EAl  2 
A21  =  EA21 
A22  -  EA22 
A31  «  EA31 
a32  «  EA32 
A41  «  EA41 
A42  «  EA42 
134$  CONTINUE 
1 349  A21  «  -A 21 
A41  «=  -a41 
GAM  «  DGAM(N) 

GAMM1  =  DGAM1  (N) 

RA  •=  DRA(N) 

RB  b  DRB(N) 

H  «  DH(N) 

C  COMPUTE  ELEMENTS  OF  E  INVERSE  FOR  THE  LAST  LAYER 
Bl 1  «  -GAM*COVA**2 
B13  «  t,/(RHO(N)*A(N)*A(N)) 

B22  «  GAMM1  *C0VA**2/RA 
B24  -  B13/RA 
B44  -  1 ,/(H*GAM) 

B33  «  -B44/RB 
B31  «  -B33*GAMM1  *H 
b42  b  1  . 

£A11-B11*A11+B1  3*A31 
EAl  2-B1 1*A1  2+B1  3*A32 
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EA21«B22*a21+B24*a41 
E A  2  2= B  2 2*A  2  2+B  2 4* a4  2 
EA  31 =B31 *Al  1+B33*A31 
EA32=B31  ,vA1  2+B33*A3? 

E  A4 1  «B4  2  *A  2 1  +  84  4*a4  1 
E  A42  -B42  *a  22+ b4  4*a4  2 

DR*  EA21  *EA32  -  EA11*EA42  -  EAl  2*EA4l  +  EA22*EA31 

D!  =  EAl 1 *EA32  +  EA21 *EA42  -  EAl  2*EA31  -  EA22*EA4l 

DENSQ  *  DR*DR  +  DI*DI 

UPMK  =  £A32*D I  -  EA42*DR 

U  PH  I  =  EA32*DR  +  EA42*D I 

V/PM  i  «  EA4l  *DR  +  EA31  *D  1 

V/PNR  =  -£A31*DR  +  EA4l*DI 

PWOTP-( ( 2,/DENSQ)*SQRTF( WPNR*WPNR+WPN I *V/PN i ) )*COVA 
PPHWP-ATAMF ( -V/PM  1  /v/PNR)- (1  ,-S  IGNF (v/PNR)  )* 

isignf(wpn; )*i  .57073 

PUDTP*  ((2./Q£NSQ)*SQRTF(UPNR*UPNR  +  UPNI*UPNi))* 

1  COVA 

PPHUP  *ATAMF(-UPMI/UPNR)-(1  .~S IGNF(UPNR) )* 
ISIGNF(UPNI)*! .57079 
PAUOV/-PWDTP/PUDTP 
dpauow=pauow/tame 

PPHUW=PPH\/P-PPHUP 
A3PH=ABSF(  PPHUV/)-3.1  41  6 
IF(A3PH)1 243,1  248,1  251 
1  251  i F ( PPHUW)1  252,1  252,1  253 
1  252  PPHUW=PPHUw+6. 2832 
GO  TO  1248 

1  253  pphuw-pphuw-6 . 2832 
1  248  I F ( PPHUW ) 395 • 395 , 396 

395  •  PPHUW-PPHUV/+3 . 1  41 6 

GO  TO  1250 

396  PPHUW-PPHUW-3.1  41 6 

1 250  PUNCH  834, PER, FREQ,GAF,PWDTP,PPHWP,PUOTPt PPHUP, PAUOW, 
1DPAUOV/.PPHUW 

834  FORMAT ( F7 . 1  ,F6.3,F7.4,5F7.2,F8.2,F7.2) 

IF  (PL0T0)740,740,745 
740  GAFL0«L0GF(GAF)/2. 302585 
IX«GAFL0*XS 
IF(PL0TR)232, 232,233 

232  PWDTP«PWDTP/(2.*SINE) 

F I Y  HLOGF ( 1 0.*PwDTP)/2, 302585 
I Y=F I Y*YS+Y01 
CALL  XYPLOT( IX, IY) 

PUDTP»PUDTP/(2.*C0SE) 

FIY«L0GF(10.*PUDTP)/2. 302585 
I Y»F I Y*YS+Y02 
CALL  XY  PL  O  f ( I X , I Y ) 

233  F I Y=L0GF(1  0.*DPAUCW)/2. 302585 
I Y«F I Y*YS+YQ3 

CALL  XYPLOTOX,  IY) 
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1X=GAFL0*X$+X02 
I F ( PLGTR) 230, 230, 221 

230  IY=*PPHW?*YPS+5.*Y04 
CALL  XYPL0T( IX, IY) 

I Y =PPHU  P*Y  PS +3 ,  *Y  04 
CALL  XYPLOT ( IX, I Y) 

231  IY=PPHUW'7PS+Y04 
CALL  XYPLOT( IX, IY) 

GO  TO  310 

72*5  IX=FREQ*XS 

IY=PWDTP*YS+Y01 
CALL  XYPLOT(IX.IY) 

IY=PUDTP*YS+Y02 
CALL  XYPL0T(IX,IY) 

IY«=PAU0W*YS+Y03 
CALL  XYPL0T( IX, I Y) 

IX=FREQ*XS+X02 
IY«PPHWP*YPS+5.*Y04 
CALL  XYPLOKIX,  IY) 

I Y  *=PPHU  P*Y  PS+3. *Y  04 
CALL  XYPLOKIX.  IY) 

IY=PPHUW*YPS+Y&4 
CALL  XYPLOT ( IX, IY) 

310  CONTINUE 

311  CONTINUE 
PUNCH  709 

709  format(//) 

C  CONTROL  CARD  NEGATIV£=NEW  ANGLE.  IF  CONTROL  CARD  =  0, 

C  A  NEW  CRUSTAL  MODEL  .  IF  CONTROL  CARD  -  A  POSITIVE 

C.  NUMBER,  THE  PROGRAM  WILL  END. 

1020  READ  7000,CNTRL 

7000  F0RMAT(1X,F9.2) 

IF(CNTRL)7001,1  030,1  021 

1021  PRINT  990 

990  F0RMAT{1X,21HEND  OF  CASE.THANK  YOU) 

END 
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GPH18  PARTOIF 


*0705 

C  THIS  PROGRAM  IS  DESIGNED  TO  CALCULATE  THE 
C  TRANSFER  FUNCTIONS  OF  THE  VERTICAL  AND  HORIZONTAL 
C  COMPONENTS  OF  LONGITUDINAL  SEIZMIC  WAVES  IN  A 
C  LAYERED  MEDIUM,  AND  THE  FIRST  PARTIAL  DERIVATIVES 
C  OF  THESE  TRANSFER  FUNCTIONS  WITH  RESPECT  TO  THE 
C  THICKNESS,  THE  ELASTIC  PARAMETERS  AND  THE  DENSITIES 
C  OF  ANY  OF  THE  LAYERS  OF  THE  SYSTEM.  BESIDES  THIS 
C  PROGRAM  CALCULATES  THE  TRANSFER  FUNCTION  FOR  THE 
C  APPARENT  ANGLE  OF  EMERGENCE,  THIS  IS,  THE  RATIO 

C  OF  THE  VERTICAL  TO  THE  HORIZONTAL  COMPONENT  AND 

c  its  first  partial  derivatives. 

DIMENSION  D(5).FL(5),A(5),RHO(5),DGAM(5),DGAM1(5), 
lDRA(5),B(5),DH(5),ELTAD(4}tELTAL(4),£LTAR(4),PAl1 
2(4),  PA12(4),PA13(4),PA14(4).PA21(4)-PA22(4). 

3PA23(4),PA24  4),PA31(4),PA32(4).PA33(4),PA34(4), 
4Pa42(4‘  Pa43(4),Pa44(4).PEA11(4),PEA12(4),PEa21(4), 
5P£A22(4),PEA31'4),PEa32(4),Pa41(4) 

D  I  MENS  I  ON  PEA4i  (4) ,  PEA.42(4)  ,DPDPR(4)  ,DPDQL  (4) , 
1DPDQR(4),DPDGL{4),DRB(5),DPDGR(4),DDRAL(4),DDRAR 
2(4) ,DDRBR(4) ,DPDPL(4),DDRBL(4) 
SIGNF(VAR)«ABSF(VAR)/VAR 

READ  CARD  TO  IDENTIFY  THE  CRUSTAL  MODEL  UNDER  STUDY 
READ  1101 
PUNCH  1101 

FORMAT (1  OH  ) 

READ  IN  THE  NUMBER  OF 
READ  2.NQL 
FORMAT? 13) 

PUNCH  2,N0L 

READ  IN  LAYER  CONSTANTS, D( 1 )«  THICKNESS  OF  THE  ITH 
LAYER.  FL(1)»PARAMETER  LANDA  OF  THE  ITH  LAYER, 
RHO(l). DENSITY  OF  THE  ITH  LAYER 
PUNCH  61 23 

6123  F0RMAT(1X,9HTHICKNESS,3X,5HLANDA,4x,7HDENSiTY,1X, 

1  10HP  VELOCITY, 2X,1  OHS  VELOCITY,/) 

DO  3  1=1  ,NOL 

READ  483, D( I ).FL( I ) ,RHO( I ) 

F0RMAT(3F10.4) 

A( I )=$QRTF(3.*FL( I )/RH0( I )) 

B(l)=SQRTF(FL(I)/RHO(l)) 

PUNCH  500,D( l),FL( I ),RHO( I ),A( I) ,B( I ) 

F0RMAT(5F10.4) 

N-NOL 


f  030 

1101 

C 


LAYERS  PLUS  ONE  OF  THE  SYSTEM 


483 


3 

500 
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C  READ  AG  I  PI e  INITIAL  ANGLE  OF  INC ID£NCE,DANG I  - 
C  INCREMENT  OF  THE  ANGLE  NQANG«NUM3ER  OF  ANGLES  TO 
C  BE  CALCULATED 

READ  4,  AG  I  PI ,DANG I,NOANG 
4  FORMAT (2F1 0.4, 13} 

C  READ  IN  INITIAL  VALl'f-'S  OF  FREQUENCY,  INCREMENT  OF 
C  FREQ,  NUMBER  OF  FREQ. 

READ  4,fmin„fin:,nof 
7001  ANG  I  P=  AG  I  PI -DANG  » 

C  READ  IN  CHANGES  IN  THE  LAYERS  PARAMETERS.  IF  THE 

C  PARTIAL  DERIVATIVE  IS  DESIRED  WITH  RESPECT  TO  THE 

c.  parameter,  hake  eltadOH.  etc.  if  not=o  no  more 

C  THAN  4  PARTIAL  DERIVATIVES  MAY  BE  CALCULATED  IN  ONE 

C  PASS  OF  THE  PROGRAM  BECAUSE  OF  MEMORY  CAPABILITIES 

C  OF  THE  1620  COMPUTER 

DO  100  |-1, HI 

read  101,  eltad(i),  eltal(i),  eltar(i) 

100  PUNCH  101,  ELTAD(I),  ELTAL(I),  eltar(j) 

101  FORMAT  (3F1  0.4) 

PUNCH  6319 

6319  F0RMAT(2X,5HFREQ.,3X,5HGAMMA,7X,2HTW,17X,2HTU,17Xf 
1  5HTW/TU,/) 

DO  310  ia«=i,noang 
FREQ=FM I N~F I NC 
5004  ANGIP-ANGIP+DANGI 

PUNCH  305, ANG I P 
305  FORMAT (F8. 3) 

C«A(N0L)/SINF(ANG IP/57.295780) 

C  COMPUTE  REUSABLE  VAR’ABLES  FOR  THE  LAYERS 
DO  1346  M«1 ,  NOL 
DH(M)  =RH0(M)*C*C 
COVA=C/A(M) 

C0V8»C/B(M) 

DRA(M)»SQRTF(A3SF(C0VA**2-1  .)) 
DRB(M)eSQRTF(A8SF(tOVB**2-1  .)) 

C0VP»1-./DRA(M) 

COVS-1 ./DRB(M) 

DPDGL(M)»2./DH(M) 

DPDGR(M)s*"2,*'FL(M)/ (DH(M)*RHO(M) ) 
DDRAL(M)»-DH(M)*C0VP/(6.*FL(M)*FL(M)) 
DDRAR(M)eC*C*C0VP/(6.*FL(M) ) 

DDRBL(M)— DH(M)*C0VS/(2.*FL(M)*FL(M)) 
DDRBR(M)«C*C*C0VS/(2.*FL(M)) 

DPDPL(M)— DH(M)*D(M)*C0VP/(6.*FL(M)*FL(M)) 
DPf'PR(M)«0(M)*C*C*C0VP/(6.*FL(M)) 
DPDQL(M)—DH(M)*D(M)*C0VS/(2.*Fl  (M)*FL(M)  ) 
DPDQR(M)=D(M)*C*C*C0VS/(2.*FL(M)) 
DGAM(M)-2.*FL(M)/DH(M) 

1 346  DGAM1  (M)-DGAM(M)-!  . 

DO  310  IFR-1,N0F 
FREQ-FREQ+FINC 
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GAF»=FREQ*D(1)/A(1) 

V/VN0=6#2831 853*FREQ/C 

Al  1“1 

Al  2=0 

A21  =0 

A22=1 

A31«0 

A32b0 

A41  =0 

A42«0 

HD«0 

COMPUTE  THE  ELEMENTS  OF  A  MATRIX. FOR  REMAINING  LAYERS 

DO  131*5  M-t,N1 

GAMoDGAM(M) 

GAMMIoDGAMI(M) 

RA«DRA(M) 

RB*=DRB(M) 

H=DH(M) 

P«4/VN0*D(M)*RA 

Q«WVNO*D(M)*RB 

SINP-SIMF(P) 

v/*»s  inp/ra 

X=RA*S IMP 
C0SP“C0SF(P) 

SINQ-SINF(Q) 

Y*»S  INQ/RB 
Z“RB*S INQ 
C0SQ»C0SF(Q) 

RH0M=RH0(M) 

ZZ-RB*C0S0 

XX=RA*C0SP 

YY«C0SQ/RA 

WW-COSP/RA 

YYY«C0SQ/RB 

DM=D(M) 

NDM-0 

B1 1  «GAM*C0SP-GAMM1  *C0SQ 
B1  2-CAMM1  *W+GAM*Z 
B1 3— (C0SP-C0SQJ/H 
Bl4»  (W+Z)/H 
B21 -  GAM*X+GAMM1  *Y 
B22— GAMM1  *C0SP-K3AM*C0SQ 
B23  — (X+Y)/H 
B24  «B1 3 

B31  «H*GAM*GAMM1*(COSP-COSO) 

B32  -H*(GAMM1  *GAMM1  *W-K5AM*GAM*Z ) 

B33  “B22 
B34  -B1  2 


153 


B4l  ^H*(GAM*AjAM*X+GAMM1  *GAMM1  *Y ) 

B42  «B31 
B43  rB21 
B44  “Bl  1 

IF(  ELTAD(M))1  0,11,10 

C  CALCULATE  THE  PARTIAL  DERIVATIVES  OF  THE  MATRIX 

C  ELEMENTS  RESPECT  TO  THE  THICKNESS  OF  THE  CORRES- 

C  PONDENT  LAYER 

1 0  ND«ND+1 
NDM=NDM+1 

PAl  1  (ND)«WVN0*(-GAM*X+GAMM1  *Z) 

PAl 2{ND)«WVN0*{GAMM1 *COS P+GAM*RB**2*C0SQ) 

Pa13(ND)«(X-Z)*WVN0/H 

PAl  4(ND) «WVN0* ( C  OS  P+RB**2  *COSQ)/H 

PA21(ND)»-WVN0*(GAM*RA**2*C0S  P+GAMM1  *COSQ) 

PA22(ND)bsWVN0*(GAMM1  *X-GAM*Z) 

Pa23(ND)=  WVN0*(RA**2*C0SP+C0SQ)  /H 
PA24(ND)-PA13(ND) 

PA  31  (NO)—  H*GAMM1  *GAM*WVNO*(X-Z ) 

PA  32  ( N  D  )  «»H*(  GAMM1  ^2*WVN0*C0SP+GAM^2*RB**2*C0SQ 
1*WVN0) 

PA33(ND)=PA22(ND) 

PA34(ND)*PA12(ND) 

PA41  (ND)=*H*(GAM**2*aA**2*VAfN0*C0SP+GAMMl  **2*WVN0 
1*C0SQ) 

PA42(ND)«PA31  (NO) 

PA^3(ND)-Pa21(ND) 

PAA1t(ND)oPAll(ND) 

11  CONTINUE 

IF(  ELTAL(M))18, 19,18 

C  CALCULATE  THE  PARTIAL  DERIVATIVES  OF  iHE  MATRIX 

C  ELEMENTS  IN  RESPECT  TO  THE  LANDA  OF  THE  RESPECTIVE 
C  LAYER 

1 8  ND«ND+l 

NDM-NDM+1 

PDPL=DPDPL(M)*WVNO 

PDQL«DPDQL(M)*WVNO 

PDGL-DPDGL(M) 

PDRAL»  DDRAL(M) 

PDRBL-  DDRBL(M) 

PAl  1  ( NO  )  =*PDGL*(  COS  P-COS  Q) -G AM*S  i NP*PD PL+GAMMl  *S  I  NQ*PDGL 
PAl  2  ( ND )  «PDGL  *(  W+Z ) +G AMM1  *(  C OS  P*PD  PL-W*PD  RAL ) /RA-KiAM* 

1  (PDRBL*S INQ+ZZ*PDQL) 

PAl  3(ND )-( S I NP*PDPL-S I NO*PDQL ) /H 
PA14(ND)^VA*/*PDPL-SINP*PDRAL/(RA*RA)+PDRBL*SINQ+ 
1ZZ*PDQL)/H 

PA21  (ND)— (X*PDGL+GAM*S  INP*PDRAL-*0AM*XX*PDPL+GAMM1  * 

1 (YYY*PDQL-(SINQ*PDRBL)/{R3*RB))+Y*PDGL) 


r>  o  o 


154 


PA22(ND)=PDGL*(C0SQ~C0SP)+PDPL*GAMM1  *$INP-GAM*SiNQ 
1*PDQL 

PA23(ND)*(SINP*PDRAL+XX*PDPL+YYY*PDQL-(SINQ*PDRBL)/ 

1  (RB*RB))/H 
PA24(ND)=PA13(ND) 

PA31  ( ND ) -H*( PDGL*( COS  P-COSQ) *( 2. *GAM-1 . )+GAM*GAMM1 * 
1(SINQ*PDQL-SINP*PDPL)) 

PA32(ND)=H*(2.*PDGL*(GAM1*W+GAM*Z)*K;AMM1  **2*(COSP* 

1  PD  PL /RA-S I N  P*PD RAL / ( RA**2 ) )  4G  AM**2*( RB*COSQ*PQQL 
2+S I NQ*PDRBL ) ) 

PA33(ND)«PA22(ND) 

PA34(ND)=Pa12(ND) 

PA41  (ND )=-+!*( 2 ,*PDGL*(GAM*X+G AMM1  *Y )+GAM**2  *( S  I NP* 

1 PDRAL+XX*PDPL )+GAMM1  **2*(YYY*PDQL-Y*PDRBL/RB) ) 
Pa42(N0)»PA31  (ND) 

PA43(ND)-PA21(ND) 

Pa44(ND)-Pa11  (ND) 

19  ’  CONTINUE 

IF(  ELTAR(M)) 22, 23,22 

CALCULATE  THE  PARTIAL  DERIVATIVES  OF  THE  MATRIX 
ELEMENTS  RESPECT  TO  THE  DENSITY  OF  THE  CORRESPONDENT 

layer 

22  ND-ND+1 
NDM-NDM+1 

?DPR«OPDPR(M)*WVNO 

PDQR-OPDQR(M)*WVNO 

PDGR=OPDGR(M) 

PDRAR«  DDRAR(M) 

PDRBR-  DDRBR(M) 

PAl 1 (ND)B,PDGR*(COSP-COSQ)-GAM*S  1NP*PDPR+GAMM1 *SINQ 
1*PDQR 

PAl  2(ND)«PDGR*W+GAMM1  *(COSP*PDPR-W*PDRAR)/RA+PDGR* 

1 Z+GAM*( POPBR*S I NQ+ZZ*PDQR) 

PAl  3(NDj  31  3/RHOM+(S INP*PDPR-S INQ*PDQR)/H 
PAl  4(ND)«-B1  4/rhom+(ww*pdpr-sinp*pdrar/(ra*ra)+pdrbr 
1*SINQ+ZZ*PDQR)/H 

1Y+GAMM1*(YYY*PDQR)-(GAMM1  *SlNQ*PDRBR/RB**2)) 

PA 21 (ND)«-(X*PDGR+GAM*S JNP*PDRAR+GAM*XX*PDPR+PDGR* 
PA22( HD ) «PDG R*  ( COS Q-COS P ) +G AMI41  *S  I N P*PDPR-G  AM* 

1 S !NQ*PDQR 

PA23(ND)=-B23/RH0M+(SINP*PDRAR+XX*PDPR+YYY*PDQR- 
1(Y*PDRBR)/RB)/H 
PA24(  ND )  =*PAl  3  ( ND ) 

PA 31  (ND)«B3’  /RHOM+H*( (COSP-COSQ)*(GAMM1+GAM)*PDGR 
1  +GAM*GAMM1 *(S INQ*PDOR~S INP*PDPR)) 

PA32(  ND )  «B  32/RHOM+H*(  2 .  *PDG  R*(GAMM1  *W+GAM*Z ) 

1 +GAMM1  **2  *(V/W*PDPR-W*PDRAR/RA)+GAM**2  * 

2(S INO*PDRBR+ZZ*PDQR) ) 

PA33 (ND)*Pa22(ND) 

Pa34(ND)«?A12(ND) 

pa4i (nD)«B41 /RH0M+H*(2,*PDGR*(GAM*X+GAMMt *Y)+GAM**2* 
1  (SINP*PDRAR+XX*PDPR)-K3AMM1**2  *(YYY*PDQR-Y 
2*PDRBR/RB)) 

PA42(ND)-PA31(ND) 
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PA43(ND)«PA21(ND) 

Pa44(ND)-Pa1 1 (MD) 

23  CONTINUE 

C  MULTIPLY  MATRICES 

£Al  1  «311*A11  +  B12*A21  +  31 3*A31  t  B1  4*a41 
EA12  -  B11*A12  +  B12*«22  +  B13*A32  +  B1 4*a42 
EA21  B21*Al 1  +  B22*A21  +  B23*A31  *  324*a41 

EA22  -  B21  *Al  2  +  B22*a22  +  B23*A32  +  B24*a42 

EA31  «  B31*A11  +  B32*A21  +  B33*a31  +  B34*a41 

£A 32  -  831*A12  +  B32*a22  +  B33*A32  +  b34*a42 

EA41  »  B41  *A1 1  +  842*A2t  +  B43*a31  +  B44*a41 

EA42  -  B41*A12  +  b42*a22  +  B43*A32  +  b44*a42 

HF-ND-NDM 
NF1 -NF+1 
IF(NDM)54,55»54 

54  DO  30  l«NF1  ,ND 

PEA1 1 ( I )«PA1 1 ( I )*A1 1 +PA1  2{ I )*A21 +PA1 3( I )*A31  + 
1 PAl 4( I )*A4l 

PEA1 2( I )=PA1 1 ( I )*A1  2+PA1  2( I )*A22+PAl  3( I )*A32+ 
1PA14(J)*a42 

PEA21  ( I )«Pa21  ( 1 )*Al  1  +PA22( I )*A21+PA23( I )*A31  + 
1  Pa24(I)*a41 

PEA22( I )«Pa21 ( 1 )*Al  2+Pa22( i )*A22+Pa23( I )*A32+ 
1PA24(I)*a42 

PEA31 ( I )=Pa31 ( I )*Al  l+PA32( i )*A21+Pa33( I )*A31  + 
1 Pa34( I )*a41 

PEA32( I )-Pa31 ( I ) *Al  2+PA32( I )*A22+Pa33( I )*a32+ 
1PA34(I)*a42 

PEA41 ( I )«Pa41 ( I )*A11+Pa42( ! )*A21+Pa43( I )*A31  + 
1Pa44(I)*a41 

30  PEA42( I )-Pa41  ( I )*Al  2+PA42( I )*A22+Pa43( I )*a32+ 

1Pa44(I)*a42 

55  CONTINUE 
All  «  EAll 
Al  2  «  EAl  2 
A21  «  EA21 
A22  «  EA22 
A31  «  ea31 
A32  «  t;A32 
A41  «  EA41 
a42  a  EA42 
IF(NF)50.1 344,50 

50  DO  31  1-1, NF 

PEAl  1  ( I  )**B1 1  *PAl  1  (  I  )+Bl  2* PA 21  ( I  )+B1  3*Pa31  ( I )+ 
1B14*PA41(I) 

PEAl 2( I )-B1 1 *PA1  2( I )+B?  2*Pa22( I )+Bl  3*PA32( I )+ 
1B14*PA42( I) 

PEA21  ( I )-B21  *PAl  1  ( I )+B22*PA21  ( I )+B23*Pa31 ( I )+ 
1 B24*PA4l ( I ) 
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PcA22( I )=62i *PA1 2( I )+B22*Pa22( i )+B2i*PA32( I )+ 

1 B24*PA42( I ) 

PEA31 ( I )-B31 *PA1 1 ( i )+B32*PA21 ( I )+B33*Pa31 ( I )+ 
1B34*Pa41(I) 

P£A32( I )=B31*PA1  2( i )+B32*Pa22( j )+B33*Pa32( 1 )+ 
1B34*Pa42(!) 

PEA41  ( I )  =341  *PAl  1  ( I  )+B42*PA21  ( I  )+b43>-?a31  ( I )+ 

1 B44*PA41 ( 1 ) 

31  PEA42( I ) ~B4l *PA1  2( I )+B42*Pa22( I )+B43*Pa32( I )+ 

1  b44*Pa42( ! ) 

'j  344  CONTINUE 

IF(ND) 60,1345,60 
60  DO  33  I =1 ,ND 

PAIKI)-PEAII(I) 

PA12(|}-PEA12MJ 
PA21 ( ! )=PEA21  ( I J 
PA22(!)=PEA22(i) 

PA31(I)-PEA31(I) 

PA32(I)«PEA32(I) 

Pa41(I)-PEa41(I) 

38  Pa42(I)=PEA42(I) 

1345  CONTINUE 
1 34?  A21  «  -A21 

a41  =  -a41 

6AH  «  DGAH(N) 

GAMM1  «  DGAH1  (N) 

RA  =  DRA(N) 

RB-DRB(N) 

H-DH(N) 

C  COMPUTE  THE  ELEMENTS  FOR  THE  E  INVERSE  OF  THE  LAST 
C  LAYER 

B1 1  «»“GAM*COVa**2 

B13  *=  1./(RHO(N)*A(N)*A(N)) 

B22  a  GAMM1  *COVA**2/RA 
B24  «  B13/RA 

B44  «  1./(H*GAM) 

B33  »  -B44/RB 

B31  a  -B33*GAMM1*H 

B42  a  1  . 

EAl  1  «B1 1  *Al  1 +B1  3*A3l 
EAl  2«u1 1  *Al  2+B1  3*A32 
EA21aS22*A21+B24^A4l 
EA22«322*a22+B24*a42 
Ea31 «B31 *A1 1 +B?3*a31 
EA32aB31 *Al  2+B33*A32 
EA41  «b42*a21 +b44*a41 
E  a42«B42*a2.2+b44*a42 

DR«  EA21*EA32  -  EA11*EA42  -  EA12*EA41  +  Ea22*Ea31 
Dl  «  EAl 1*EA32  +  EA21*EA42  -  EA12*EA31  -  EA22*EA4l 
DENSQ  «  DR*DR  +  Dl*Dl 
UPNR  «  EA32*DI  -  EA42*CR 
UPN I  «  EA32*DR  +  EA42*DI 

UDTP  a  ((2./DENSQ)*SQRTF(UPNR*UPNR  +  UPN l*UPN l ))*COVA 
PHUPD«ATANF(-UPNI/UPNR)-(1  .-S IGNF(UPNR) )* 


ISIGNF(UPNI)*  1  .57079  -  ' 

V/PNI  =  EA4l*DR  +  EA 31  *0  I 
V/PNR  =  -EA31*DR  +  EA4l*DI 

V/DTP  »((2./DENSQ)*SQRTF(WPNR*V/PNR+V/PNI*V/PNI )  )*C0VA 
AU0VH7DTP/UDTP 

PHWPD --A7ANF ( — WPN  I  /V/PNR )“( 1  . -S  IGNF(WPNR)  )* 

1 S  IGNF ( V/PN  I  )*1  .57079 
PHUOV/=?:-iV/PD-PHUPD 
ABPH-ABSF ( PHU0W)-3. 1  41  6 
IF(ABPH)6l 0,61 0,611 

611  IF(PHU0v/)6l  2.61  0,6l  3 

612  PHUOV/ -PHUOW+o.  2832 
GO  TO  610 

613  PHU  OW-PHU  0V/-6 . 28  3 2 

610  PUNCH  303.  FPvEQ#GAF#WDTP,PHV/PD,UDTP,PHUPD#AUOW,PH'JOW 
303  FORMAT ( 2Fo .4, 6F9.4) 

C  FOR  INVERSE  MATRIX  OF  LAST  LAVER 
DO  32  I «1 # MD 
PA21  ( I )«-Pa21 ( I ) 

PA41 ( I )=-Pa41 ( I ) 

PEAl  l  ( I )=Bl  1  *PAl 1 ( I )+31  3*PA31 ( I ) 

PEAl  2( I )=B1 1  *PA1  2( I )+Bl  3*PA32( I ) 

PEA21  ( j )=B22*Pa21 ( I )+B24*Pa41  ( I ) 

PEA22( I )=B22*Pa22( I )+B24*Pa42( I ) 

PEA31  ( I )~B31 *PAl  1  ( I )+B33*Pa3j ( I ) 

PEA 32 ( I )  =»B31  *PAl  2(  I  )+B33*Pa32(  I ) 

PEA41  ( I J-B42*Pa21 ( I )+B44*PA4l  ( I ) 

PEa42( I )=b42*Pa22( I )+B44*Pa42( I ) 

PDR=>PEA21  ( I  )*PEA32(  I  )-PEAl  1  ( I  )*PEA42(  I  )- 
1  PEAl  2(  I  )*PEA41  (  I  )+PEA?-2(  I  )*PEA31  ( I ) 

PD  I  *=PEAl  1  ( I  )*PEA32(  I  )+PEA21  ( I  )*PEA42(  I  )- 
1  PEAl  2( I )*PEA31 ( I )“PEA22( I )*PEA4l ( I ) 

PUR=€R*PEA42( I )+D l*PEA32( I )-EA42*PDR-EA32*PD I 

PUI=D!*PEA42( I )-DR*PEA32( I )-EA42*PD l+EA32*PDR 

PWPv=~DR*PEA31  ( I  )-D  l*PEA4l  ( I  )+EA4l  *PD  I+EA31  *PDR 

PWI «DR*PEA41 ( I )-D l*PEA31 ( I )-EA4l *PDR+EA31  *PD I 

D2R=DR*DR-DI*DI 

D2l=»2.*DR*Dl 

PU  PNR=PU  R*D 2R+PU I  *D  2 1 

PU  PN I “PU I *D  2R-PU  R*D 2 1 

PD  ENSQ=D  2R*D 2R+D  2  i  *D  2 1 

PUDTP“(2./PDENSO)*COVA*SQRTF(PUPNR*PUPNR+PUPNI*PUPNI ) 
PPHUPD“ATANF(-PUPNI/PUPNR)-(1  .-SIGNF(PUPNR))* 

1 S  IGNF(PUPN  I  )*1  .57079 
PWPNR=PWR*D 2R+PW 1  *D  2 1 
PV/PN I  =~PV/R*D  2 1 +PW I  *D  2  R 

PV/DTP=(2./PDENSQ)*C0VA*SQRTF(PWPNR*PV/PNR+PWPNI*PWPNI) 

PPHV/PD=ATANF(-Pv/PNI/PWPNR)~(1  .-signf(pwpnr))* 

1 S IGMF(PWPN I )*1  .57079 

PRR«-F.a42*PEa31  ( I  )+EA32*PEa41  ( I  )+PEA42(  I  )*EA31- 
1EA4l*PEA32(l) 

PRI«Ea32*PEA31 ( I 5+EA42*PEA4l ( I )-EA4l*PEA42( I )- 
1 EA31 *PEA32( I ) 

P2R«Ea42**2-EA  32**2 
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P2!«2.*Ea42*EA32 
PRPNR»PRR*P2R“PR I *P2 I 
PRPNI«PRR*P2i+PRI*P2R 
PRDESQ«P2R*P2R+P2 I *P2 1 
PMUOW-0 ,/PRDESQ)*SQRTF( PRPNR-PRPNR+PRPN l*PRPN i ) 
PPHUOW=ATANF(-PRPNi/PRPNR)-U IGNF{ PRPNR) )* 
1SIGNF(PRPNI)*1 .57079 

PUNCH  307, I ,FREQ,GAF, PWDTP, PPHWPD, PUDTP, PPHUPD, PMUOw 
1.PPHUOW 

307  FORHAT( 12,2F8.4,6F9.4) 

32  CONTINUE 
310  CONTINUE 

C  CONTROL  CARD  TO  CONTROL  RECYCLE 

1020  READ  7000,CNTRL 
7000  F0RMAT(F10.2) 

IF(CNTRL)7001 ,1030,1  021 

1021  CONTINUE 
END 
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APPENDIX  lit 


**  FERNANDEZ  GPH  18  DlAGR 

C  THIS  PROGRAM  1$  DESIGNED  TO  PLOT  THE  PARTICLE 

C  MOTION  DIAGRAMS  OF  THE  GROUND  MOTION  CORRESPONDING 

C  TO  AN  INCIDENT  LONGITUIDINAL  WAVE  THROUGH  A 
C  LAYERED  SYSTEM,  AT  DIFFERENT  VALUES  OF  THE 

C  PARAMETER  GAMMA 


DIMENSION  PWDTP(l6),PUDTP(l6),PHWOU(l6) 

14  PUNCH  TAPE  100 

100  FORMAT  (1H«) 

IX*  3800 
IY*  0 

CALL  XYPLOT( IX, IY) 

1X*0 

IY«2500 

CALL  XYPL0T(!X,IY) 

DO  3  LG*1  ,1  o 

C  READ  THE  VERTICAL  AMPL  I TUDE( PWDTP)  THE  HORIZONTAL' 

C  AMPLITUDE  (PUDTP)  AND  THE  PHASE  D I FF ERENC E ( PHWOU ) 

READ  2  , PWDTP (LG) , PUDTP (LG), PHWOU (LG) 

2  F0RMAT(20X,F7.2„7X,F7.2,22X,F7.2) 

IF  ( PWDTP(LG) )4,1  0,4 

4  CONTINUE 

DO  20  L*1 ,3 
READ  21,  ZRO 
21  FORMAT (20X,F7. 2) 

20  CONTINUE 

3  CONTINUE 

1 0  NE*LG-1 

LG-0 

DO  7  JY *4,20,8 
G, =2400- JY*1  00 
DO  6  JX*5,33,7 
0X«JX*1  00 
LG*LG+1 
DO  200  NS*1 ,7 
FNS=NS 

IX*0X-400.+FNS*1  00. 

IY*OY 

200  CALL  XYPLOT(IX,IY) 

DO  201  MS-1,7 

!X«OX 

FMS-MS 

iy«qy-4oo.+fms*i  00. 
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201  CALL  XvPL0T(!X,IY) 

IF  (NE  .G)12,11,11 

11  00  5  i  **! ,  1 2 
F I « I  -1 
W-FI/2. 

IX»PUDTP{LG)*100.*C0SF(W+PHW0U(LG)+3.1415)+  OX 
|  Y  =— PV/DT  P  ( LG )  *1  00,*C0SF(W)  +  OY 
[»  CALL  XYPL0T(IX,1Y) 

3  CONTINUE 

7  CONTINUE 

12  CONTINUE 

READ  500,  MORE 
500  FORMAT( 12) 

IF  (M0RE)13, 13,14 

13  PRINT  15 

15  FGRMAT(21H£ND  OF  RUN,  THANK  YOU) 

END 
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APPENDIX  IV 


SEISMIC  ANALYZER  PROGRAM 


MAIN  PROGRAM 

The  main  program 

1)  Reads  the  scale  factors  of  the  different 
graphs  of  the  output. 

2)  Reads  the  identification  of  the  earthquake  and 
the  characteristics  of  the  spectral  analysis  to 
be  performed,  such  as  data  window  length,  digi¬ 
tization  interval,  frequency  resolution  of  the 
spectrum,  number  of  passes  of  the  record  through 
the  digital  filter  and  the  limits  of  the  fre¬ 
quency  of  the  spectrum. 

This  part  of  the  program  controls  the  use  of  ail  the 
subprograms,  which  are  optional  and  can  be  omitted  If  de¬ 
sired. 

Finally,  this  main  program  reads  the  frequency  ampli¬ 
tude  and  phase  response  of  t'ne  three  instrumentr  of  the 
station  and  the  digitized  values  of  the  seismograms  of  the 
three  components. 

The  program  can  be  recycled  to  perform  the  operations 
with  another  set  of  records . 

FILTER  SUBPROGRAM 


This  is  an  optional  subprogram  whose  use  is  controlled 


by  the  Main  Program. 
Its  object  is: 
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1)  To  remove  the  zero-shift  of  the  record  by  a  cor¬ 
rection  based  on  the  average  of  all  the  data 
points . 

2)  To  remove  the  linear  trend  that  could  be  intro¬ 
duced  in  the  records  by  mislocation  of  the  zero 
line.  In  practice  this  step  traces  a  straight 
line  through  the  points  and  relocates  all  the 
points  about  this  line. 

3)  To  apply  the  linear  digital  filter  indicated  be¬ 
fore  . 

WINDOW  SUBPROGRAM 

This  is  an  optional  subprogram  whose  use  is  controlled 
by  the  main  program. 

Its  object  is: 

1)  To  weight  the  seismic  record  in  such  a  way  that 
amplitudes  rear  the  first  arrival  of  the  P  wave 
are  emphasized  and  .later  arrivals  minimized.  This 
is  done  according  to  the  "hamming"  operation  de¬ 
scribed  before.  The  time  interval  of  the  window 
is  controlled  by  the  main  program  and  is  optional 
for  each  case. 

2 )  To  smooth  out  the  spectrum. 

3)  To  plot  the  seismogram  after  these  operations 


have  been  performed 


At  this  point  the  seismogram  of  the  P  wave  is  ready 
for  Fourie:  integration. 

FOURIER  INTEGRAL  SUBPROGRAM 

This  is  an  optional  subprogram  that  can  be  used  in  two 
different  forms:  a)  cosine  Fourier  integral  for  power  spec¬ 
trum  calculation  from  autocorrelation  data,  or  b)  complex 
Fourier  integral.  This  optional  control  is  given  in  the 
Main  Program. 

Its  object  is: 

1)  To  calculate  the  Fourier  integral  of  the  records. 
The  numerical  Integration  is  performed  using 
Fllon's  method. 

2)  To  correct  the  amplitudes  and  phases  of  the 
Fourier  analysis  by  the  frequency  response  of 
the  Instrument  for  amplitudes  and  phase. 

3)  To  print  an  output  of  the  individual  spectra  of 
the  three  components .  Though  this  output  Is  not 
necessary  for  the  object  of  this  dissertation, 
nevertheless  it  can  be  used  in  related  investiga¬ 
tions  . 


RATIO  SUBPROGRAM 

This  is  an  optional  program  whose  use  is  controlled 
by  the  Main  Program. 

Its  object  is: 

1)  To  divide  the  vertical  component  of  ground  motion 
(as  represented  by  the  output  of  the  previous 


program)  by  the  horizontal  component.  The  result 
is  printed  and  plotted  on  a  logarithmic  scale. 

2)  To  find  the  phase  difference  between  the  vertical 
and  the  horizontal  components  at  the  different 
frequencies  of  the  Fourier  analysis.  The  results 
are  printed  and  plotted  on  a  logarithmic  scale. 

In  addition  to  the  above,  the  program  gives  the  total 
vector  amplitude  of  ground  motion  at  the  frequencies  of 
the  Fourier  analysis. 

By  a  convenient  control  the  program  may  be  recycled 
to  analyze  several  records  in  succession. 


Ib5 


**  FERNANDEZ  GPH  18  SEISMIC  ANALYZER. MAIN 

C  THIS  PROGRAM  IS  D£S iGNED  TO  DIVIDE  THE  SMOOTHED 

C  SPECTRA  Of  THE  VERTICAL  AND  HORIZONTAL  COMPONENTS  OF 

C  THE  P  V/AVES  AND  TO  PLOT  THE  TANGENT  OF  THE  APPARENT 

C  ANGLE  OF  EMERGENCE  VERSUS  FREQUENCY. 

C  THE  PROGRAM  IS  COMPOSED  OF  A  MAIN  PROGRAM  AND  FOUR 
C  OPTIONAL  SUBPROGRAMS. 

DIMENSION  Y(  900)tFM0D(3,  60) ,GA I N( 3) ,COR(1 00) , 

1  XRS  ( 3}  •  XRO(  3) , YRS  (3)  ,YRO(  3)  #XFS(  3) . XFu(  3) ,  YFS  ( 3) 
2,YFO(35,XCS(3),XCO(3),YCS(3),YCO(3),PHA(3,  60) 

COMMON  Y , FMOO  , NT , COR, FM 1 N , F I NC , NOF ,G A i N , H, FNFL , 

1 PLOTF  ,XCS,XCO,YCO#SMOTH#M#TM,YCS,XFS,XFC,YFS,YFO, 
2XSS,XS0,YSS,YS0,NAC,TL,C0MP, JC, INDS,PHA,T1 ,T2,T3, 
3T4fT5,T6,CT1,CT2,CT3,CT4,CT5,C76,TF1 .TF2,TF3,TF4, 
4TF5,TF6,CTF1 ,CTF2,CTF3,CTF4,CTF5,CTF6 
C  READ  THE  SCALE  FACTORS(XS, ETC)  AND  ORIGIN  POINT 
C  (01 , ETC)  FOR  THE  PLOT. 

DO  270  1=1 ,3 

READ  1  ,XRS( !),XRO( I ),YRS( I ),YR0( I ).XFS( I ),XFO( f ) 

READ  1  ,YFS(I),YF0(I),XCS(I),XC0(!),YCS(I),YC0(!) 

1  FORMAT (oFI  0.0) 

270  CONTINUE 

READ  411 ,XSS,XSO#YSSfYSO 
411  F0RMAT(4F10.0) 

C  '  READ  A  CARD  WITH  THE  OPERATIONS  AND  SUBPROGRAMS 

C  TO  BE  USED. 

C  IF  FNFL=>  POSITIVE  NUMBER  THE  FILTER  PASS  IS  USED. 

C  IF  PLOTR£=  POSITIVE  NUMBER,  THE  DIGITIZED  RECORD  IS 

C  PLOTTED.  IF  PLOTF I =  POSITIVE  NUM3ER  THE  FILTERED 

C  RECORD  IS  PLOTTED.  IF  SMOOTH^  POSITIVE  NUMBER  THE 
C  SUBPROGRAM  WINDOW  IS  USED.  IF  FOURT=POS IT ! VE  NUMBER, 

C  THE  SUBPROGRAM  FTRAN  IS  USED,  IF  V/URAT=POS  IT  I  VE 
C  NUMBER,  THE  SUBPROGRAM  RATIO  IS  USED, 

READ  31.  FNFL,PLOTR,PLOTF,SMOTH,FOURT,WURAT 
31  FORMAT(bFIO.O) 

C  READ  A  CARD  WITH  THE  PARAMETERS  OF  THE  SPECTRUM 
C  COMP*  POSITIVE  NUMBER  FOR  FOURIER  COMPLEX  INTEGRAL 
C  FMIN=MIN I MUM  FREQUENCY, F  C«INCREMENT  OF  THE  FREQUENCY 

C  ,NOF=NUMB£R  OF  POINTS  TO  BE  CALCULATED, NAC*POS IT  I VE 
C  NUMBER  FOR  PLOT  OF  SMOOTH  RECORD. 

C  INDS=POSTI VE  FOR  LISTING  OF  INDIVIDUAL  SPECTRA. 

READ  30,FMIN.FINC,N0F,NAC, INDS,COMP 
30  F0RMAT(2F1  0.3,313, FI  0.0) 

PUNCH  TAPE  3 
3  F0RMAT(1H») 
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C 

22 

2 

C 

C 

C 

C 

C 

C 

4 

C 

C 

C 

C 


283 


99 

284 


C 

C 

24 

6 

7 


9 


10 

8 

12 

13 

16 

15 

18 

17 

5 

20 


IDENTIFICATION  OF  SPECTRUM 
READ  2 


PUNCH  2 
F0RMAT(49H 

GAIN  OF  RECORDS  AT  CENTRAL  FREQUENCY, 1=Z,2=N-S, 

3«E“W  COMPONENTS.  TL  LENGTH  OF  THE  WINDOW  IN  SECONDS 
CALC=0  IF  NOT  CALIBRATION  CURVE  IS  USED  TO  FIND 
GROUND  MOT  I  ON, CALC =1  IF  THE  RESPONSE  CURVE  IS  USED 
THE  RESPONSE  CURVE  OF  THE  INSTRUMENTS  IS  GIVEN 
WITH  THE  DATA 


READ  4.GA1N0  5,GAlN(2),GAlN(3),K1  ,TL,CALC 
F0RMAT(3F10.0,F1  0.7,2F10.3) 

READ  THE  MAGNIFICATION  ,CT1,ETC,THE  PHASE  DELAY, 
CTF!,ETC  AND  THE  CORRESPONDING  P£RI0DS,T1,  ETC', 

FOR  SIX  SELECTED  POINTS  OF  THE  RESPONSE  CURVE. 

INTERMEDIATE  VALUES  OF  THE  PERIOD  ARE  INTERPOLATED 

DO  5  JC=1,3 

IF  (CALC)  283,284,283 

READ  99,Ti,CT1,T2,CT2,T3,CT3 

READ  99,T4#CT4,T5,CT5,t6,CT6 

READ  99,TF1 ,CTF1 ,TF2,CTF2,TF3,CTF3 

READ  99. TF4,CTF4, TF5 ,CTF5 , TF6, CTF6 

format(6fi  0.2) 


CONTINUE 

DO  6  1=1,1600/4 

READ  24, Y ( I  ),Y(  1+1  ),Y(  l+2).Y(  1+3), Y(  1+4)  Y(  1+5) 
1Y(l+6),Y(l+7),Y(l+8),Y(l+9>,Y(l+10),Y(l+l1),Y(l+12 
2),Wl+f3) 

THE  LAST  CARD  OF  EACH  COMPONENT  MUST  BE  A  CARD 

WITH  9999 

F0RMAT(14F5.0) 

IF(Y(l)-9999.)6,7,6 

CONTINUE 


NT=l-2 


h~hi 

IF(PLOTR)  8,8,9 
DO  10  1=1, NT 
F  1  =  1—1 

I X=F I *H*XRS ( JC )+XRO( JC ) 
i Y *=*Y ( I  )*YR$i  UC)+YRQ(  JC; 
CALL  XYPLOT(IX,IY) 
CONTINUE 
IF(FNFL)1 3.1 3,1 2 
CALL  FILTE(FNFL) 

IF  (SMOTH)1 5,1 5,1  6 
CALL  WINDO(COR) 

IF  (FOURT)  17,17,18 
CALL  FTRAN  (COR) 
CONTINUE 
CONTINUE 

JF(WURAT)  19,19,20 
CALL  RATIO(FQ) 
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19  PUNCH  TAPE  3 

READ  34,  MORE 
34  FORMAT(FIO.O) 

IF  (MORE)  21,21,22 
21  PRINT  23 

PUNCH  23 

23  F0RMAT(20HEND  OF  RUN, THANK  YOU) 

PAUSE 
END 


FERNANDEZ 


GPH  18  FILTER  SUB. 


C 

C 

C 


463 


464 

470 


462 


465 


SUBROUTINE  F  ILTE(FZ) 

THIS  IS  AN  OPTIONAL  SUBROUTINE  CONTROLED  BY  THE 
MAIN  PROGRAM.  ITS  OBJECT  IS  TO  DETREND  THE  RECORD 
AND  TO  FILTER  IT. 

DIMENSION  Y(  900),FM0D(3,  60) ,GA I N(3) ,COR( 1 00) 
lXRS(3)#XR0(3),YRS(3)fYR0(3),XFS(3).XF0(3),YF$(3) 
2,YFO(35,XCS(3),XCO(3),YCS(3),YCO(3),PHA(3,  60) 
COMMON  Y,FMOD,NT,COR,FMIN,FINC,NOF,GAlN,H,FNFL, 
1PL0TF  ,XCS,XCu,YCO,SMOTH ,M,TM,YCS*XFS,XFO. YFS, YFO, 
2XSS,XSO.YSS,YSO,NAC,TL,c6MP,JC, INDS,PHA,T1  ,T2,T3, 
3T4,T5,T6,CT1,CT2.CT3,CT4,CT5#CT6,TF1.TF2,TF3,TF4, 
4TF5,TF6,CTF1 ,CTF  2,CTF3,CTF4,CTF5*CTFo 
NH=NT/2 


FNL=NH 

YL=0. 

DO  463  I  =1  ,NH 
YL*=YL+Y(  I ) 
YL=YL/FNL 
YR«0. 

NJ-NH+1 


ENR-NT-NH 
DO  464  I  «*NJ,NT 
YR-YR+Y( I ) 

YR==YR/FNR 

TO-C. 

NTR-1 

NFL«FNFL 

I F(NFL)461 ,461,462 

FNT-NT 

NF=NT-5 

C1=2.*(YR-YL)/FNT 

C2«(YR-YL)/2. 

N=0 

DO  465  l-1,NF,2 

F I  «l 

N°N+1 

Y(N)«(Y ( I  )+Y ( l+5)+2;*(Y ( 1+1  )*>Y(  i+4))+3,* 
1  ( Y (  I+2)+Y(I+3)))/12.-(FI+2.5)>C1+C2 
NT«N- 
FT*2.5*H 
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TC.-TO+FI 

H=2t*H 
Y  ( N+1  }  =0. 

Y(N+2)=»0. 

IF(NFL-NTR)466,466,467 

467  NL=N/2 
YL  =0 

DO  463  |«1,NL 

468  YL  -YL+Y  ( i ) 

NR1<=NL+1 

YR=0, 

DO  469  |«NR1,N 

469  Y  R=Y  R+Y ( 1 ) 

FNL«=NL 

FNR=N-NL 

YL=YL/FNL 

yr=yr/fnr 

NTR=NTR+1 
GO  TO  470 
466  FT-TO 

IF(PLOTF  )115, 115,116 
116  DO  471  1=1, N 
F I  =  I  “1 

I X=(F l*H+FT)*XFS( JC)+XFO( JC) 
IY=Y( i )*YFS( JC)+YFO( JC) 

CALL  XYPLOT(!X,!Y) 

471  FT=FT+H 
461  CONTINUE 
115  CONTINUE 
RETURN 
END 


**  FERNANDEZ  GPH18  WINDOW  SUB. 

C  THIS  IS  AN  OPTIONAL  SUBROUTINE  CONTROLED  BY  THE 

c  main  program,  its  object  is  to  pass  the  record 

C  THROUGH  A  DATA  WINDOW  OF  VARIABLE  LENGTH. 

SUBROUTINE  WJNDO(HAM) 

DIMENSION  Y(  900).FM0D(3,  60),GAIN(3),C0R(1 00), 
1XRS(3).XP0(3),YRS(3),YR0(3)#XFS(3),XF0(3),YFS(3) 
2,YFO(33,XCS(3),XCO(3),YCS(3),YCO(35,PHA(3,  60) 
common  y,fmod.nt,cor,fmin,finc,nof,gain,h,fnfl, 

1 PL0TF  ,XCS,XCO,YC0,SM0TH,M,TM,YCS,XFS,XF0,YFS,YFQ, 
2XSS,XS0,YSS,YS0#NAC,TL,C0MP,JC, INDS.PHA.T1  ,T2,T3, 
3T4,T5,T6,CT1,CT2tCT3,CT4,CT5,CT6,TF1 .TF2,TF3,TF4, 
4TF5,TF6#CTF1  ,CTF2,CTF3,CTF4,CTF5,CTF6 
C  TOTAL  LENGTH  OF  THE  REC0RD=TL 
DO  43  1=1, NT 
BI=M 
T=H*B1 

IF  (T-TL)  83,83,84 


16Q 

83  i.  =0.54+0.4G*C0SF(3.l4l5S*T/TL) 

GO  TO  85 

84  WD=0. 

85  Y(I)=Y(1)*WD 

C  IF  PLOT  OF  THE  RECORD  IS  DESIRED  LET  NAC=1 

IF  (NAC)  81,81,82 
82  iX=T*XCS(JC)+XCO(JC) 

I Y  =*Y (  I  )*YCS(JC)+YCO(JC) 

CALL  XYPLOT( IX, I Y } 

81  CONTINUE 

43  CONTINUE 

RETURN 
END 

**  FERNANDEZ  GPH  18 

SUBROUTINE  FTRAN(R) 

C  THIS  IS  AN  OPTIONAL  SUBROUTINE. 

C  FOURIER  ANALYZE  THE  RECORDS  AND  TO  CONVERT  TO 

C  GROUND  MOTION  THE  AMPLITUDE. 

DIMENSION  Y  (  900)tFM0D(3,  60)fGAlN(3)  ,COR(1  00) , 
1XRS(3),XR0(3),YRS(3),YR0(3),XFS(3).XF0(3),YFS(3) 
2,YF0(3),XCS(3),XC0(3),YCS(3),YC0(3),PHA(3,  60) 
COMMON  Y,FMOD,NT,COR,FMIN,FINC,NOF,GAlN,H,FNFL, 

1 PLOTF  , XCS,XCO,YCO,SMOTH,M,TM, YCS,XFS, XFO,YFS, YFQ, 
2XSS,XS0tYSS,YS0,NAC,TL.C0MP,JC, INDS.PHA.T1  ,T2,T3, 
3T4,T5,To,CTl ,CT2,CT3,CT4,CT5,CT6,TF1 ,TF2,TF3,TF4, 
4TF5 , TF6, CTF1 , CTF2, CTF3, CTF4, CTF5  #  CTFo 
S I GNF ( VAR) =ABSF ( VAR) /V AR 
IF(INDS)1  002,1  002,1  003 
1003  PUNCH  1  001,  JC 

1001  FORMAT (26h  SPECTRUM  FOR  COMPONENT, 2X, II ,//) 
PUNCH  1000 

1000  FORMAT ( 4 3H  FREQ.  PERIOD  AMPLITUDE 

1  PHASE  ) 

1002  IDE=(NT-2)/2 
IDE-2* IDE 
FREQ=FMIN-F INC 
DO  3  1=1, NOF 
FREQ=FREQ+F INC 
PER=1  ./FREQ 
W=6.283l8*FREQ 
SWDT=SINF(W*H) 

CWDT=COSF(W*H) 

ZETA«=W*H 

ALFA«(ZETA**2+ZETA*SWDT*CWDT-2.*(SWDT**2))/ZETA**3 
BETA=2.*(ZETA*( 1 ,+CWDT**2)-2.*SWDT*CWDT ) /ZETA**3 
GAMA *4 . *( SWDT-Z  ETA*CWDT ) /Z  ZT A**3 
SWT=0. 

CWT-1 . 
s=0. 

C“Y(1  )/2. 

DO  1  J«"2,  IDE, 2 


FOURIER  INTEGRAL 
ITS  OBJECT  IS  TO 


---  xr 
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N=J 

SA=SWT 

SWT=SWT*CWDT+SWDT*CWT 
CWT«CV/T*CWDT-SA*SWDT 
c=c+gama*y( J)*CWT 
IF(C0MP)4,4,5 

5  S=S+GAMA*Y(J)*SWT 

4  CONTINUE 

$A*SWT 

S1'  r  — SV/T  *CV/D  T +S  WD  T  *CWT 

CWT«CV/T*CWDT-SA*SWDT 

IF(C0MP)6,6,7 

7  s=s+beta*y( J+1  )*SWT 

6  CONTINUE 

1  C=C+BETA*VY(  J+1  )*CWT 

SA=SWT 

SV/T  =3  WT  *CWD  T +S  WD  T*CWT 
C WT=C WT *C WD  T-S  A*S  WD  T 
C=C+Y( IDE+1  )*CWT/2. 

S-S+Y ( IDE+1  )*SWT/2. 

IF(C0MP)8,8,9 

8  FM0D(JC,I)«2.*H*(Y(IDE+1  )*$WT*ALFA+C) 

FMOD ( JC , I ) =SQRTF ( FHOD ( JC# I ) ) 

PHA(JC,l)=0. 

GO  TO  10 

9  C=Y ( IDE+1  )*SWT*ALFA+C 
S=~(Y( IDE+1  )*CWT+Y(1  ) )*ALFA+S 

FHOD( JC, I )«SQRTF( (S*S+C*C)*#1 591549)*H 
PHA(JCf l)=ATAN(-S/C)-(1  ,-S IGNF(C) )*S IGNF(S)*1  .5707963 
C  CORRECT  FOR  INSTRUMENT  RESPONSE 

I F(CALC) 1 1 9,1  0,1 1  9 

119  T=PER 

IF(T2-T)1  20,1  20,1  21 

121  CTI=CT1+(CT2-CT1 }*(T-T1 )/(T2-T1 ) 

GO  TO  140 

120  IF  (T3-T)1  22,1  22,1  23 

123  CT I =CT2+( CT3-CT2) *(T-T2 ) /(T3-T 2) 

GO  TO  140 

127  IF(T4~T)1 24,1  24,1  25 

1 25  CTI=CT3+(CT4-CT3)*(T-T3)/(T4-T3) 

GO  TO  140 

124  iF(T5-T)1 26,1  26,1  27 

1 2?  CTI®CT4+(CT5-CT4)*(T-T4)/(T5-T4) 

GO  TO  140 

1  26  CTI=CT5+(CT6-CT5 )*(T-T5 )/(T6-T5) 

140  CONTINUE 

FMOD ( JC , I ) «FMOD ( JC , I ) *1  000 . /CT I 
IF  (T2-T)150, 150,151 

151  CTF !=C7F1  +(CTF2-CTF1  )*(T-TF1  )/(TF2-TF1 ) 

GO  TO  160 

150  IF(T3*"T)1 52,1 52,1 53 

153  CTFI«CTF2+(CTF3-CTF2)*(T-,F2)/(TF3-TF2) 

GO  TO  160 
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152 

155 

154 

157 

156 
1 60 

C 

1C 

170 

15 

1 6 
13 
3 


IF  (l4-T)1 54f 1 54,1  55 
CTFI-CTF3+(CTF4-CTF3)*(T-TF3)/(TF4-TF3) 
GO  TO  1  60 

IF  (T5-T)156, 156,157 

CTF I =CTF4+( CTF5-CTF4 ) *( T-TF4 ) / (TF5~TF4 ) 
GO  TO  160 

CTF  I  =CTF5+( CTF6-CTF5 ) *( T-TF^ ) / ( TF6-TF5 ) 
PHA( JC, I )=PHA(JC, I }“CTF I 
GO  TO  170 

NORMALIZE  VALUES  TO  VERTICAL  COMPONENT 


FMOD(JC, i )=FMQD(JC, I ) *GA I N( JC) /GA I N(1 ) 

CONTINUE 

I F ( INDS )1 3,1 3,1 5 

PUNCH  l6,FREQ.PER,FM0D(JC,l),'>HA(JC,l) 
FORMAT (4F1  2.4; 

CONTINUE 

CONTINUt 

RETURN 

END 


** 
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SUBROUTINE  RAT IO(FQ) 

C  THIS  SUBROUTINE  HAS  AS  ITS  OBJECT  TO  DIVIDE  THE 

C  SPECTRUM  OF  THE  VERTICAL  COMPONENT  BY  THE  VECTORIAL 

C  SUM  OF  THE  SPECTRA  OF  THE  HORIZONTAL  COMPONENTS. 
DIMENSION  Y(  yOO),FMOD(3,  60),GAIN(3),C0R(1  00), 

1  XRS( 3) ,XRO( 3) , YRS( 3) ,YRO( 3) ,XFS ( 3) ,XFC( 3) ,YFS( 3) 
2fYFO(3),XCS(3),XCO(3),YCS(3),YCO(3),PHA(3,  60) 

COMMON  Y , FMOD , NT, COR, FM I N, F I NC, NOF.GA I N , H, FNFL , 
1PL0TF  ,XCS,XCO,YCO,SMOTH,M,TM,YCS,XFS,XFO,YFS,YFO, 
2X$S,XS0,YSS,YS0,NAC,TL,C0MP,JC, INDS,PHA,T1 ,T2,T3, 
3T4,T5,Tb,CT1  ,CT2,CT3,CT4,CT5,CT6,TF1  .TF2,TF3,TF4, 
4TF5,TF6,CTF1,CTF2,CTF3,CTF4,CTF5,CTFb 
S I GNF ( VAR) -=ABSF ( VAR ) /VAR 
PUNCH  6324 

6324  FORMAT (  2X,43HSPECTRUM  OF  THE  APPARENT  ANGLE 
1  OF  EMERGENCE,/) 

PUNCH  6325 

6325  F0RMAT(2X,6HPEPI0D,5X,5HFREQ.,4X,7HM0DULUS,5X, 
15HPHASE,2X,11H.0TAL  ampl.,/) 

DO  1  I  =1  ,NQF 
R  J=  I  “1 

FREQ=FMIN+RJ*F INC 
PER=>1  /FREQ 

HM0D=SQRTF(FM0D(2,  ; )**2+FMCD(3,  l)**2) 

W0U=FMGD(1 ,  I  )/HMOD 

PR«FM0D(2, I )*C0SF(PHA(2, I ) )+FM0D(3, I )*C0SF(PHA(3, I ) ) 
PI -FMOD (2. I)*SINF(PHA(2, l))+FM0D(3, I )*S INF(PHA(3, I)) 
PHAH-ATAN(-PI/PR)-(1  .-SIGNF(PR))*$IGNF(PI)*1  .5707963 
DPHA-PHAO  ,  I  )-PHAH 
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A3PH=ABSF(DPHA)-3.1416 

if(abph)6io,6io,6ii 

611  !F(DPHA)6l2#6lO,6l3 

612  DPHA=0PHA+6.2832 
GO  TO  610 

613  DPHA=DPHA-6.2832 

61  0  AMPT  =  SQRTF ( 7  MOD (1,1) **2+HM0D**2 ) 

PUNCH  2,  PER, FREQ, V/OU ,DPHA. -MPT 
2  FORMAT (4f1 0.4, FI  0.2) 

FREQ=L0Gr( FREQ )/2. 302585 
IX=FREQ*XSS+XSO- 
WOU=LOGF(WOU)/2. 302585 
iY=V/OU*YSS+YSC 
CALL  XYPLOTOX,  IY) 

1Y=DPHA*1  00. +2 100. 

CALL  XYPL0T( IX, IY) 

1  CONTINUE 

RETURN 
END 
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